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ABSTRACT 

We examine the physical basis for algorithms to replace mixing-length theory (MLT) in stellar evolu¬ 
tionary computations. Our 321D procedure is based on numerical solutions of the Navier-Stokes equa¬ 
tions. These implicit large eddy simulations (ILES) are three-dimensional (3D), time-dependent, and 
turbulent, including the Kolmogorov cascade. We use the Reynolds-averaged Navier-Stokes (RANS) 
formulation to make concise the 3D simulation data, and use the 3D simulations to give closure for the 
RANS equations. We further analyze this data set with a simple analytical model, which is non-local 
and time-dependent, and which contains both MLT and the Lorenz convective roll as particular sub¬ 
sets of solutions. A characteristic length (the damping length) again emerges in the simulations; it is 
determined by an observed balance between (1) the large-scale driving, and (2) small-scale damping. 

The nature of mixing and convective boundaries is analyzed, including dynamic, thermal and com¬ 
positional effects, and compared to a simple model. We find that (1) braking regions (boundary 
layers in which mixing occurs) automatically appear beyond the edges of convection as defined by the 
Schwarzschild criterion, (2) dynamic (non-local) terms imply a non-zero turbulent kinetic energy flux 
(unlike MLT), (3) the effects of composition gradients on flow can be comparable to thermal effects, 
and (4) convective boundaries in neutrino-cooled stages differ in nature from those in photon-cooled 
stages (different Peclet numbers). The algorithms are based upon ILES solutions to the Navier-Stokes 
equations, so that, unlike MLT, they do not require any calibration to astronomical systems in order 
to predict stellar properties. Implications for solar abundances, helioseismology, asteroseismology, 
nucleosynthesis yields, supernova progenitors and core collapse are indicated. 

Subject headings: stars: evolution, oscillations, supernovae; convection; turbulence 


1. INTRODUCTION 

Make everything as simple as possible , but no simpler. 
-Albert EinsteinO 

Stars contain three dimensional (3D), turbulent 
plasma. They are much more complex than the sim¬ 
plified one dimensional (ID) models we use for stellar 
evolution. Computer power is not adequate^ at present 
for well-resolved (i.e., turbulent) 3D simulations of whole 
stars for evolutionary timescales. 

We attempt to tame this complexity by (1) use 
of 3D simulations as a foundation, (2) application of 
the R eynolds-Averaged Navier - Stokes (RANS) pr oce- 
dure dMeakin fc Arnett! I2007bl ; IViallet, et al.l 1201311 to 
these simulations to discover dominant terms (closing the 
RANS system), and (3) construction of simple physical 
models, consistent with the 3D simulations, for use in 
stellar evolution codes. We call this approach “321D” 
because a central feature is the projection of 3D simu¬ 
lations down to ID f or use as a replaceme nt for mixing- 
length theory 1MLT: [Bohm-Vitcnsdll95~8lh The process 
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is designed to allow testing, extension, and systematic 
improvement. 

Formally, the RANS equations are incomplete unless 
taken to infinite ordei0: they must be closed by trunca¬ 
tion at low order to be useful. This need for truncation is 
due to the nature of the Reynolds averaging, which allows 
all fluctuations rather than only dynamically consistent 
ones. Closure requires additional information to remove 
these extraneous solutions. Using 3D simulations avoids 
this problem by providing only dynamically consistent 
fluctuations. 

As a complement to the full RANS approach, we con¬ 
sider approximations which focus on dynamics; these 
provide a connection to historical work on convection 
in astrophysics and meteorology. Such a minimalist step 
may be easier to implement in stellar evolutionary codes, 
and still provide physical insight. In the turbulent cas¬ 
cade, kinetic energy and momentum are concentrated in 
the largest eddies. Our approximate model contains both 
the largest eddies and the Kolmogorov cascade. 

1.1. Historical Background 

Erika Bohm-Vitense developed the version of mixing- 
length theory used in stellar evol ution in the 1950s 
(jVitensel 119531 iBohm-Vitensd 1195811 . prior to the pub¬ 
lication in the west of A ndrey Ko l mogor ov’s theory of 
the turbulent cascade (lKolmogorovlll962lf . MLT might 

8 This occurs because the momentum equation is nonlinear, so 
that eac h level of correlation requir es the next higher level for its 
solution (Tritton 1988; Vallis 2006), giving an infinite regression. 
See also ICubarsil (120101) . 


























2 


Arnett et al. 


have been different had she been aware of the original 
work (lKolmogorovHl94lT) . Edward Lorenz showed that 
a simple c onvective r oll had chaotic behavior (a strange 
attractor. iLorcnd 119631 ). Ludwig Prandt l developed the 
theory of boundary layers (jPrandtl fcTiet]ens[jl934(l . as 
well as the original version of MLT (lPrandtlfll925f l. All 
these ideas will be relevant to our discussion, which is 
based, as far as possible, upon experimentally verified 
turbulence theory and 3D simulations, and free of astro¬ 
nomical calibration. 
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Figure 1 . The 3D turbulen t energy casc ade on a logarithmic scale 
of sizes; see lDavidsonl d 20041) ; IPopef (12000) . The arrow indicates the 
direction of net energy flow. The range of applicability of Implicit 
Large Eddy Simulations (ILES) and Direct Numerical Simulation 
(DNS) are shown. See text for discussion. 


The 3D turbulent energy cascade is illustrated in Fig¬ 
ure [T| The turbulent motion is driven at the largest 
scale (the “integral” scale), which contains most of the 
kinetic energy. These motions are unstable and break 
up into smaller-scale flow patterns dominated by inertial 
forces (the “inertial subrange”). This continues to scales 
small enough for microscopic effects (viscosity) to finally 
provide damping of the flow at the Kolmogorov scale. 
Both the inertial subrange and the dissipation range are 
insensitive to the details of the boundary conditions at 
the integral scale, and are “universal” in this sense. We 
use the term “universality” to mean the property of in- 
sensitivity to boun dary conditions at the integral scale. 
iKolmogorov ( 194 11) found the striking result that the rate 
of dissipation is insensitive to the value of the viscos¬ 
ity, but is determined by the rate that the largest-scale 
flows feed the cascade. This behavior of the non-linear 
flow “hides” the microscopic value of the viscosity. We 
use Kolmogorov theory to describe the flow in the range 
where universality holds. 

Direct Numerical Simulations (DNS) resolve the small 
scales at which dissipation happens, and can extend up 
to the inertial range, but not to stellar scales. Implicit 
Large Eddy Simulations (ILES) can extend from stellar 
(integral) scales down to the inertial range, but not to 
the dissipation range. Fig. [l] illustrates both. 

Landau objected to the notion of complete universal¬ 
ity on the grounds that the largest scales were subject 
to boundary con ditions which would be specific to th e 
case in question (lLandau fc Lifshiti 119591 IFrischl 119951) . 
We will incorporate this idea by splitting the turbulent 
flow into two parts: the integral-scale motion and the 


turbulent cascade. As an aid to understanding the in¬ 
tegrated properties of the integral-scale motion, we are 
guided by the simplest model of a convective roll, due to 
iLorenzI (119631 ). This model contains the famous Lorenz 
strange attractor, and exhibits chaotic behavior. It also 
agrees surprisingly well with three-dimensional (3D) sim¬ 
ulations of turbulent convectio n associated with oxygen 
burning prior to core co llapse (|Meakin fe Arnettll2007bl ; 

I Arnett fe Meakinll2011bll . This approximation does lack 
multi-mode behavior, as compared with the simulations, 
w hich are dominated by fiv e low order modes (see Fig. 1 
in lArnett fc Meakinll2011bl ): this may affect the accuracy 
of the representation of intermittency at large scales and 
of coherent structures. 

Our challenge is to simplify this very complex problem, 
with time dependence and an astronomically large num¬ 
ber of degrees of freedom, down to a feasible level for use 
in a stellar evolutionary code, without losing important 
features. Our approximation, 321D, is an attempt to in¬ 
crease physical realism at feasible cost in computational 
complexity. It is desirable to avoid astronomical calibra¬ 
tion as far as possible, and base changes upon behavior 
quantified in laboratory and numerical experiments. In 
particular, we do not validate our approximation by how 
well it reproduces standard MLT results. By basing ap¬ 
proximations on 3D ILES simulations that (1) exhibit 
turbulence, (2) have non-uniform composition, and (3) 
resolve dynamic boundary behavior, it is possible to re¬ 
move some of the vagueness inherent in many theoretical 
treatments of convection. 

We will compare the global properties of turbulent con¬ 
vection from numerical and analytical viewpoints in Sec¬ 
tion^ examine the structure and nature of boundaries of 
convection zones in Section [3l and summarize our conclu¬ 
sions in SectionjdJ In an appendix we provide a derivation 
from 3D fluid flow equations for some useful expressions. 

2. GLOBAL BEHAVIOR OF CONVECTION 

lArnettl _ (119941): iBazan fe Arnettl (1199811 : 

lAsida fc Arnett] ([200(1( 1 found that 2D simulations 
of stellar oxygen burning developed large fluctua- 
tions at the boundaries of t he co nvective region. 
iKuhlen. Wooslev. fe Glatzmaierl (12003H found that 3D 
simulations of the same stage g a ve no such boundary 
fluctuations. iMeakin fc Arnettl (12006H did both 2D 
and 3D simulations and showed that the discrepancy 
was due to a d ifferent choice of boun dary condition: 
IKuhlen. Wooslev. fc Glatzmaicd (12003H used rigid 
boundaries at the edge of the convective region, while 
the other simulations included dynamically-active stable 
layers surrounding the convection, a more realistic 
choice. Nevertheless, all obtained a convective velocity 
of u ~ 10'cm/s. The global character of the velocity field 
seemed to be insensitive to the details of the convective 
boundary, although these fluctuations are an important 
part of the physics of the boundary itself (and the extent 
of the convective region). This insensitivity allows us to 
separate the g lob al pro blem from the boundary problem 
(see also ICanutol 119921) : in this section we focus on the 
global problem. 

The turbulent kinetic energy equation may be inte¬ 
grated over a convective region; in the steady state limit 
this gives a global balance between driving on the inte¬ 
gral scale, and dissipation at the Kolmogorov scale (see 
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Fig. HD. This balance has been verified experimentally 
and numerica ll y as a common feature of tu rbulence (e.g., 
iTennekes fe Lumlevl Il972t iDavidsonl I2004T) . This intro¬ 
duces a length scale, the depth of the convective zone, 
into the problem. 

2 . 1 . The Turbulent Cascade 

Using _a_c lassica l radiative viscosity 

(iMihalas fc Mihalasl [1984 1. the Reynolds number is 
Re ~ 10 8 at the base of the solar convection zonc0- 
Numerical simulations and laboratory experiments 
become turbulent for roughly Re > 10 I * 3 , so fluid flows 
in stars are strongly turbulent if, as we assume for the 
moment, rotational and magnetic field effects may be 
neglected. 

For homogeneous, isotropic, and steady -state turbu¬ 
lence, the Kolmogorov relation (jFrischll 1995T) between the 
dissipation rate of turbulent kinetic energy per unit mass 
et, velocity i>, and length scale £ is 

U = (1) 

I Arnett. Meakin. k, Young! (12009( 1 found that e* = 
0.85 v 3 ms /£ cz , where £ cz is the depth of the convective 
zone, and v rms is the average convective velocity; see 
their Eq. 6 and nearby discussion, and references to other 
studies which report such coefficients. Fo r homogeneous, 
isotropic turbulence, iKolmogorovl (1194111 predicted a co¬ 
efficient 4/5 for a region well away from boundaries. This 
factor of 0.8 might change for a strongly stratified region, 
which would have flow better described by plumes than 
convective rolls. 

Eq.[l]is a global constraint, averaged over fluctuations, 
and applies to each length scale A in the turbulent cas¬ 
cade, so 

e t ~ (Av\) 3 /X, (2) 

for all scales A, or, 

Av\ ~ (e t A)T (3) 

so that the velocity variation across a scale A is Av\, 
which increases as As. The largest scales have the largest 
velocities, and are dominated by advective transport 
(macroscopic mixing). 

The velocity gradient across the scale A is 

Av\/X ~ ef /\i, (4) 

and increases with decreasing A. The smallest scales have 
the largest velocity gradients, and are eventually domi¬ 
nated by microscopic mixing (ionic diffusion, radiative 
diffusion, and viscosity). A description of the cascade 
needs both large and small scales; Eq. [3] implies that the 
largest (integral) scales have most of the kinetic energy 
and momentum, while Eq. [4] implies that the smallest 
scales have the fastest relax ation times, w hich is con¬ 
sisten t with simulations (e.g.. lArnctt. Meakin. fc Young! 
2009). 


2.2. Limitations of Resolution 

lLandau fe Lifshitzl (U959I) . §32, estimated the number 
of degrees of freedom in a region of turbulent flow to 
be A ~ (Re) 9 / 4 . Laminar flows with free boundaries 
become unstable at roughly Re ~ 10 3 * . A direct numer¬ 
ical simulation (DNS) would require well over 10 8 zones 
to resolve the cascade for this marginally unstable case. 
Using Re ~ 10 8 (see Section [241 . implies a need for more 
than 10 18 zones for the Sun, far beyond current computer 
capacity. 

There may be a smarter way. Kolmogorov’s great in¬ 
sight is that turbulence hides the details of the viscous 
dissipation by the nonlinear interactions of the cascade, 
so that the dissipation rate is determined by macroscopic 
parameters. Simulations show a multimode behavior 
dMeakin fc Arnett! 12007 b). but only N ~ 5 dominant 
modetH for ~ 10 8 zones. This dramatic reduction in 
complexity suggests the use of implici t larg e eddy sim¬ 
ulations (ILES, see Fig. |T| and iBorisI 12007 ) which ap¬ 
proximate small scale behavior by a Kolmogorov cas¬ 
cade. Our approach is to assume that this simplification 
holds for very large Reynolds numbers, and to exam¬ 
ine the consequences. Simulations which are presently 
feasible have effective Reynolds numbers limited by nu¬ 
merical resolution, but are sufficiently high to give truly 
turbulent solutions. State of the art simulations, with 
both improved algorithms a nd more powerfu l co mput¬ 
ers, support this approach dPorter fc Woodward! 120001 : 
iHerwig. et al.ll20l4 iCampbell. et al. 1120151) . 


2.3. Dynamics: MLT to 321D 

As an aid to the reader, Table Q] gives the correspon¬ 
dence of selected variables in three different theoretical 
approaches to turbulent convection: MLT, the Lorenz 
model, and the RANS formulation. MLT is ID (radial), 
the Lorenz model is 2D (radial and transverse), while the 
RANS analysis is 3D projected to ID. MLT is static, the 
Lorenz model and the RANS equations are time depen¬ 
dent. MLT is local (no spatial derivatives of velocity) 
while the Lorenz model is mildly nonlocal (it uses global 
derivatives over the roll), and the RANS equations are 
non-local. Comparison of MLT and Lorenz gives a sense 
of transverse versus radial properties. 

In MLT the buoyant acceleration is approximately in¬ 
tegrated over a mixi ng length £mlt to obtain an aver¬ 
age velocity u (e.g.. IVitense 119531 iBohm-VitenselH958I : 
iKippenhahn fc Weigertlll99C 1. 


u 2 =gf3 T A v(|^f). (5) 

The superadiabatic excess AV is defined in Table Q] and 
112.41 Here g is the gravitational acceleration, /3t = 
— (<9 In p/<9 In T)p is a thermodynamic variable (for uni¬ 
form composition; see ^2.41 for the nonuniform case), Hp 
is the local pressure scale height, and £mlt is an ad¬ 
justable length scale (the mixing length). 

Eq. [5] requires that AV > 0 for the velocity u to be 
a real number. The velocity depends only on the local 


9 Using only a classical plasma viscosity due to ion collisions, the 
Reyn olds number would be even larger <1Arnett, Meakin &; Vial let! 
[20TD . 


10 See lHolmes, Lumlev, &; Berkooj (119961) and more recent work 
on principle component analysis and other techniques which at¬ 
tempt to exploit the reduction in complexity. 
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Table 1 

Correspondence of some variables in MLT, Lorenz and RANS 


quantity 0 

MLT 6 

Lorenz 6 

RANS d 

comment 

dissipation length 

(■ a 2 /8)H P 

— 

id ~ 0.8 lez 

icz is convection zone depth 

Meakin & Arnett (2007b) 

horizontal 0 gradient 

A V = V-V e 

(^)T 3 /To 

— 

Smith & Arnett (2014) 

radial gradient 

V e — V a 

(^)T 2 /T 0 

— 


imposed gradient 

Vr Va 

(2Rp)Ti/T 0 

— 


convective velocity 

Eq. [5] 

U 

u' 

algebraic (MLT) versus ode 
local (MLT) versus nonlocal 

turbulent heating 

none 

ignored or 
u 2 \u\/e d 

{{n'f)/t d 

Arnett, Meakin, & Youne (2009) 

kinetic energy flux 

assumed 
cancellation 
by symmetry 

assumed 
cancellation 
by symmetry 

(p'u'u ■ u/2) 
no cancellation, 
asymmetry 

Meakin & Arnett (2007b, 2010) 

buoyancy flux 

u/3j-^AV 

I/3 t9 «T 3 /To 

o 

1 

MLT ignores composition gradients 

enthalpy flux 

puCpTAV 

lpC P uT 3 

pCp(u'T') 

Viallet, et al. (2013) 

acoustic energy flux 

none 

none 

{P’u’) 

small for low-mach flow 

composition 6 flux 

undefined 

none 

p(Y’u') 

Arnett (T996) 

Y e flux 

undefined 

none 

P(Y/W) 



a The MLT variables are all defined in the radial direction. RANS projects a 3D average onto the radial direction. 
The Lorenz model has b oth radial and horizontal gradients (lArnett fc M cakin 2011b; Smit h fc Arnetdl2014D . 
b Smith_^Arnetj (|201^). 

c Arnett, Mcakin, & Young (20091: lArnett fc Meakinl (I2011bf) : t is the roll diameter. 
d Meakin fc Arnettl ll2007bll : Viallet. et al.l (I2013f) . 
e lArnetd (I1996H . Y = Y e + E«Yi = Y e + 1/A and Y e = YiZiYi. 


value of the superadiabatic gradient AV. There are ob¬ 
vious problems with regions in which such integration 
extends past a boundary. 

There have been a num be r of att empts to gen¬ 
eraliz e MLT; e.g., lUnnol (119611) . IGough (119671 
I1977T). lArnettl (I1969H . iStellingwerfl (fl976l). Kuhfusd 


Xiond ( 1986 1. iDeng. Bressan &; Chiosi ( 1996 1 


Xiong, et al. 1997 1. _ Hansen fe Kawalerl ( 1994 1 


Deng, Xiong, & Chan (120061 1. etc. Working backward, 
Eq. [5] may be expressed as a co-moving acceleration 
equation for a vector field u: 


du/dt = B — V, (6) 

where B is a gener alized driving term and T > a corre¬ 
sponding drag term dPrandtl fc Tietieiisl (119341) , Ch. V). 
A hydrostatic background will be assumed; see Appendix 
i jAl Similar equations result from (1) study of the nonlin¬ 
ear development of the Rayleigh-Taylor instability (RTI), 
and from (2) applications of Reynolds-Averaged Navier- 
Stokes (RANS) analysis to 3D simulations of turbulent 
convection. 

If the driving is due to buoyancy alone, (see H2.4I for 
nonuniform composition), —g(Sp/p) ~ g/3 T AV, then 
B « g/SyAV. If the drag is represented by V s=s u/r, 
where t = l d /\u\, then we have 

du/dt = du/dt + (u • V)u = g/3 T AV — u/r. (7) 


This is basically a statement of Newtonian me chanics, 
with driving by buoyancy and damping by drag. iGoughl 


(1977) gives a historical c ontext going back to iPrandtll 
( 1925 1 and to iBiermaiml (|1932ll . The early attempts, 
and many of the recent ones, have used a kinetic theory 
model, in which the mixing length was a sort of mean free 
path. In contrast, we interpret Eq. [6] as a model of the 
momentum equation for fluid dynamics, involving struc¬ 
tures such as waves, convective rolls, or plumes. Because 
it is non-local, Eq. [7] allows formally stable regions to be 
convective, unlike MLT, because of finite velocities. This 
may be relevant for composition mixing in weakly stable 
regions, and the mass contained in convective regions. 

Taking the dot product of Eq. [7] with u gives a kinetic 
energy equation, 

d(u 2 /2)/dt = u • g/lrAV — u 2 /r, (8) 

for which the steady-state solutiorQ is Eq. [5j with t d = 
i 2 MLT /&Hp, and AV > 0. In Eq. [51 negative values 
of AV are allowed; t h is permits buoyant deceleration 
(iBrummell. Glune. fc Toomrdl2002l ). The singula rities in 
MLT at the convective zone bound aries (§9 in IGough 
II9771) . and in boundary layers (§40 in lLandau & Lifshitzl 
11959ft are removed 1 ^!. 

11 Care must be taken (for negative u ) with the sign of the transit 
time t and the deceleration. 

12 The singularities in this case occur in Prandtl’s equations for 
a boundary layer as the velocity perpendicular to the surface goes 
to zero. In a star the motion does not go to zero but becomes 
wave-like rather than turbulent. 
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The flow is relative to the grid of the background stel¬ 
lar evolution model, so the co-moving time derivative of 
turbulent kinetic energy leads to 

d(u 2 /2)/dt = c?t(u • u)/2 + V • F K , (9) 

where Fk = pu(u • u)/2 is a flux of kinetic energy. The 
generation of the divergence of a kinetic energy flux in 
this way is robust for dynamic models; it occurs in the 
more precise RANS approach (Eq. US] as well as Eq. [SJl. 

We may write Eq. [U] as 

d t {u ■ u)/2 + V ■ F k = u {B-V). (10) 

In a steady state, the divergence of turbulent kinetic en¬ 
ergy flux is zero only if there is a local balance between 
the driving and the drag terms. Otherwise turbulent ki¬ 
netic energy flux may be non-negligible. The turbulent 
kinetic energy flux smooths the distribution of turbulent 
kinetic energy between regions in which it is generated in 
excess, and the whole turbulent region. The drag term 
is usually relatively smooth in comparison to the driv¬ 
ing term, which can be strongly peaked. Turbulent ki¬ 
netic energy transport is especially important if convec¬ 
tion is driven by cooling near the photosphere, so that the 
(negati ve) buoyancy is localized and the stratification is 
strong. iMeakin fe Arnett) (12010|1 have shown that strat¬ 
ification enhances the asymmetry in convective kinetic 
energy flux for driving from t he top, and reduces it fo r 
driving at the bo t tom; s ee also lStein fc Nordlundl (|1989ft : 
ICattaneo. et al. I ()1991l ). This asymmetry is small for 
shallow convective zones, growing with stratification. 

This behavior does not occur in MLT, which en¬ 
forces an exact symmetry between up-flows and down¬ 
flows so that V ■ Fk = 0. Although simulations 
of 3D atmospheres exhibit strong downward (neg¬ 
ative) net fluxes of kinetic energy, such informa¬ 
tion was not i nclude d in MLT fits for such atmo¬ 
spheres (iTramnedachl 120071 : iTrampedach fc Steinl 120111 : 
iMagic. et al.ll2014f). Simulation s of 3D red-giant atmo- 
spheres bv lLudwig fc Kucinskad (120121 1 indicate that the 
fits to MLT require at least a two parameter family, as 
have sim ulations of deeper c onvection. In the red giant 
model in lViallet. et all ( 2013 1. the downward directed ki¬ 
netic enei^y^u2LT§ackf s 35% of the maximum enthalpy 
flux. iStcin fc Nordlundl (119981) find that their solar model 
has a downward directed kinetic energy flux which is 10% 
of the enthalpy flux. This downward kinetic energy flux 
must be compensated for by a larger (outward) enthalpy 
flux. This kinetic energy flux is accompanied by a mo¬ 
mentum flux, which affects the convective boundary, as 
shown in d3.81 These are nontrivial differences relative 
to MLT, and may have implications which are detectable 
with asteroseismology as deviations from the predictions 
of MLT models. 

At present, stellar evolution theory has no turbulent 
heating term. This is inconsistent with Kolmogorov 
theory, which states that turbulent kinetic energy is fed 
back into the thermal bath at the rate given by Eq. [Q 
From the viewpoint of a dynamic model (e.g., Eq. [6]), 
this is a “frictional” cost of moving e nergy by convec¬ 
tion. lArnett. Meakin. fc Younal ( 2009 ) show that ener- 

13 Alternatively one might take the view that this is included 
in the MLT “convective flux” by construction, but this conflates 
different physical effects. 


getic self-consistency requires that the usual stellar evo¬ 
lution equations must be modified to include such a heat¬ 
ing term, or equivalently, to explicitly include terms for 
heating by buoyancy work an d divergence of ki n etic en¬ 
ergy and acoust ic fluxes (see Arnett. Meakin. fc Younel 
120091 Ea. 20-22: iMocak. et al.ll20l4 821.5. S21.6). The 
Kolmogorov term appears as heating in the internal en¬ 
ergy equation and cooling (damping) in the turbulent 
kinetic energy (acceleration) equation. Total energy is 
conserved; turbulent kinetic energy is transformed into 
heat. 

It may be more convenient to apply the heating term 
directly, rather than use the buoyancy work and diver¬ 
gence of turbulent kinetic energy and acoustic fluxes, as 
the velocity is available from solution of Eq. [TJ Turbu¬ 
lent heating (and divergence of kinetic energy flux) may 
have implications for the standard solar model and so¬ 
lar abundance^. Such heating may also be important 
for the motion of convective burning shells into electron- 
degenerate fuel. 

In the local, steady-state, limiting case, the left-hand 
side of Eq. [8] vanishes, and an equation similar to 
Eq. [5] results, but with a turbulent damping length in¬ 
stead of a mixing length. In simulations this is the 
lesser of the de pth of the convective zone or 4 pres¬ 
sure scale height (j Arnett fc Meakinll2011bD . With this 
chang e, the cubic equation o f Bohm- Vitense may be de¬ 
rived (ISmith fc Arnettll2014D . and we recover a form of 
MLT. 

Had it been available, Bohm-Vitense might have iden¬ 
tified the mixing length with the Kolmogorov damping 
length (Eq. [TJ) . However, Kolmogorov found the damp¬ 
ing length £d to be the depth of the turbulent region, so 
that it is not a free parameter, unlike MLT. There is a 
further issue: e t is the average dissipation rate, not the 
instantaneous local value (u 3 /ld) which fluctuates ove r 
time and space (see Fig. 4 in IMeakin <fc Arnettll2007bl) : 
that is, u ^ v except on average. This is reminiscent of 
the RANS approach (1 12.61 and 112.71) . 

Suppose we assume that the integral scale motion is 
that of a 2D convective roll, where du/dt is given by 
Eq.0 Using this and a corresponding thermal energy 
equation, we obtain a form of the classic Lorenz equa¬ 
tions, but with a nonline ar damping term provid ed by 
the Kolmogorov cascade (lArnett fc Meakinll2011bD . Be¬ 
cause of the time lag, as implied by the time needed to 
traverse the cascade from integral to Kolmogorov scales, 
the modified equations are even more unstable than the 
original ones, and have chaotic behavior\^\ 

14 lArnetti d20l?) suggested that the flux of turbulent kinetic 
energy was simply responsible for a change in radiative luminosity 
in the solar model. The situation is more complex. The finite 
negative luminosity of turbulent kinetic energy flow is compensated 
by an increased positive enthalpy flux, and a radiative flux. This 
modifies the therm al str ucture. The turbulent momentum flux in 
the braking region (833 extends the well-mixed region beyond the 
conventional Schwarzschild estimate; these effects would modify 
the solar model in the same sense. 

15 This upper limit to the turbulent damping length may be 
related to increasing stratification. The development of plumes 
and their Rayleigh-Taylor instability will enhance the turbulent 
drag, reducing the increase in t&\ see 823 

16 Direct integration shows that, even for no time lag in dissi¬ 
pation, chaos sets in slowly at a Reynolds number Re between 600 
and 700. 
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2.4. Nonuniform Composition 

In Eq. [7] it was assumed that the density fluctua¬ 
tion which drives the buoyancy could be represented by 
— g(Sp/p ) ~ g/3 T AV, involving only a fluctuation in tem¬ 
perature. This is only true for uniform composition and 
mild stratification. The formulation makes use of the 
expansion of pressure fluctuation, 


p' = (dP/dT) p . Y T' ± (, dP/dp) T , Y p' 


+{dP/dY) TtP Y', 

(11) 

which may be written as 


P'/P = -Pt{T’/T) - Py(Y'/Y) 


+(P/ps 2 )P'/P , 

(12) 

where 


Pt = ~ (d In p/d In T)y,p, 

(13) 

Py = —(dlnp/d\nY)T,p, 

(14) 

s 2 = (, dP/dp) T ,Y■ 

(15) 


Here s is the sound speed. The composition vari- 
able Y denote s the number of free particles per baryon 
(jArnettl I1996D . and is e ssenti ally the in verse o f the 
mean molecular weight u , (|Kippenhahn fc \\eicertlH990l : 
lHansen fe Kawalerlll994f) . An illustrative and simple ex¬ 
ample is the ideal gas, P = IZpYT. For subsonic flows, 
\P'/P\ ~ ( u/s ) 2 , where u/s is the Mach number of the 
flow, and is srnalfH. In MLT, the pressure fluctuation is 
assumed zero (no acceleration by pressure dilatation), so 

Pt(T'/T) ± Py{Y'/Y) « -(p'/p), (16) 

and it is further assumed that Y' = 0 to obtain Eq. 0 
Even in the limit of negligible pressure fluctuations, vari¬ 
ations in Y enter in a way similar to variations in 
T, so even small composition variations can be signifi¬ 
cant when superadiabatic temperature variations are also 
small. Many of the difficulties found using MLT are re¬ 
lated to situations in which Y' ^ 0: overshooting, semi¬ 
convection, and entrainment. 


2.5. Dynamics: Rayleigh-Taylor instabilities 

There seems to be a deep connection between Eq. [71 
Rayleigh-Taylor instabilities (RTI), and turbul ent mix¬ 
ing. An almost identical equation (Eq. 4.1 in lAbarzhil 
l201Gh is used to describe the nonlinear development 
of the RTI into the turbulent mixing regime. Unlike 
canonical Kolmogorv turbulence, the RT turbulent 
mixing is statistically unsteady, and involves the 
transport of potential and kinetic energies as well as 
enthalpy. Because of its importance in a variety of high 


ener 

2002 

?y-aensiiy (rirtu j 
; Kane, et al.J 

conditions (Zjeiaovicn naizer 
1997k iRemington. et all 1999k 

Dimonte, et al.| |2004 

: Remington. Drake. & Rvutovl 


imental effort for its study as well as an extensive 
literature have developed. 

The RTI happens when a heavier fluid overlays a 
lighter one, proceeding from linear instability of pertur- 


17 Near boundaries the approximation P'/P ~ 0 fails because 
pressure fluctuations p rovid e the transverse acceleration necessary 
to divert the flow; see ES 


bations (lChandrasekharlfl961lh to mildly nonlinear mo¬ 
tion of bubbl es and spikes, and then to nonlinear turbu¬ 


lent mixing (lAbarzhil 12010T) . The initial acceleration is 


one-dimensional, but as instability develops, the motion 
breaks symmetry and approaches isotropy (as seen in a 
co-movi ng frame), much like the cascade in steady tur¬ 
bulence (I Frischl 119951) . The essential difference between 
stellar convection and RTI is that the RTI is not con¬ 
tained, while convection operates within a definite and 
slowly varying volume. This means that the vertical and 
the transverse scales are causally conne cted in convec¬ 
tion, but may be independent in the RTI (lAbarz MM). 

Inconsistency between experimental and numerical in¬ 
vestigation of the RTI in the nonl inear regime led to the 
on, problem (|Dimonte. et al.ir2004D . The RTI in the limit 
of strong mode-coupling can be initiated to have self¬ 
similar evolution, so that the amplitude (diameter of the 
bubble Db (x hb) evolves as hb ~ a^Agt 2 , where A is the 
Atwood number (density ratio. [Chandrasekhar! 1 196~H) . g 
is gravity and t the elapsed time. The simulation value 
a.b ~ 0.025±0.003 is smaller than the experimental value 
ab ~ 0.057 ±0.008. This discrepancy seems to have been 
resolved by the idea that unquantified errors in the exper¬ 
imental initial conditions were the cause. To the extent 
that such uncertainties cannot be precisely known, this 
suggests a statistical approach, and illustrates the need 
for combined theoretical, experimental, and numerical 
st udies. _ 

iMeakin fe Arnettl (l2007bl ) found that regions of their 
simulated convection zone exhibited recurring “bursts” 
of convection (see their Fig. 4). These bursts, although 
multi-modal (n ~ 5), see m to share the chaotic behav¬ 
ior of the iLorend (Il963fl model of a single-mode con¬ 


vective roll ( Arnett fc Meakinll2011bD . This encourages 


the use of Eq. [71 wh ich is related to the momentum- 
driven model of RTI ( lAbarzhil l2010f) . for timescales less 
than or of order of the transit time. For longer, evolu¬ 
tionary timescales (stellar convection) we need to aver¬ 
age over fluctuations, which means averaging over sev¬ 
eral transit times for the convective roll (see Eq. [TB] be¬ 
low). These bursts result from underlying physics sim¬ 
ilar to that in the RTI; their short timescale behavior 


T-ITIC 

stoc] 

chanism. Arnett & Meakiri 2(1 

lH or equivalently, 

iastic excitation of oscillations. 

Goldreich & Kumar 

199C 

; IGoldreich. Murray. & Kumar 

1994 Aerts. et al. 


Mai . 


2.6. Filtering Fluctuations 

The weak coupling between driving at the large scale, 
and dissipation at the small scale, allows time dependent 
fluctuations of significant amplitude in luminosity and 
turbulent velocity. The term du/dt (Eq. [6] and Eo. fT71) 
is needed for chaotic fluctuations and wave generation. 
These fluctuations have a cellular structure in space and 
time; if there are many cells, with random phases, the 
fluctuations in t he average total luminos ity are reduced 
by cancellation (I Arnett fc Mcakinll2011bD . 

Fluctuations are fundamental features of turbulence 
and mixing. Because of sensitivity to initial condi¬ 
tions which can never be known with complete accuracy, 
descriptions of turbulence should be statistical in na¬ 
ture, even though the equations are deterministic dFrischl 
119951) . Turbulent simulations can be said to be numer- 
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ically converged only in a statistical sense. Eventually 
trajectories will diverge. Lyapanov exponents character- 
ize this divergence , a feature characteristic of turbulence 
(|Mannevilld l20ld ) which makes turbulent mixing so ef¬ 
fective. Unlike the diffusion picture, in which a stellar 
mixing front moves radially, limited by the random walk 
of mean-free-path strides, turbulent mixing involves a 
network of trajectories throughout the space of the tur¬ 
bulent region, laced with inhomogeneities, which finally 
disappear at the Kolmogorov scale. 



6 , , R 7 , 

radius / (10 cm) 



Figure 2. Fluctuations in velocity vg (in cm/s) versus radius. 
The top panel shows a sequence of snapshots, one for each time- 
step of St ~ 0.5s. The lower panel shows the running average over 
300 such steps (150 s), starting at 12 times, each separated by 20 
steps (10 s). The vertical scales are identical. Averaging gives 
linear cancellation of small time scale fl uctuations, but long er time 
scale variations survive. See text and IViallet. et all (20131 ). model 
OB. The simulation shows a flow involving several (~ 5) prominent 
modes which decay and reform. See text for details. 


In stratified regions, mass conservation constrains the 
flow, but it tends to change the cross-sectional area of the 
plumes as opposed to limiting their range. Although the 
flow is locally wild with fluctuations, these tend to can¬ 
cel upon horizontal and time averaging, leaving a much 
more placid behavior due to the cancellation of random 
phases. Fig. [2] illustrates this for a particular but repre¬ 
sentative case; the velocity in the theta direction, vg, is 
shown as a fu nction of radius, f rom the oxygen burning 
data set in IViallet. et al.l (|2013l ). The top panel shows 
the instantaneous value of vg (in units of cm/s) for a se¬ 


quence of time steps St ~ 0.5s. The bottom panel shows 
the running average (a horizontal average, i.e., over a 
spherical surface of radius r) of the same variable over 
300 such time steps (150s), stepping forward over 20 time 
steps (10s) at a stride, on the same velocity scale. The 
amplitude in the bottom panel is much reduced by can¬ 
cellation; what does remain is the larger length scale, 
as suggested by the cascade idea discussed in ■ nrn The 
cancellation does not work for quadratic terms; they re¬ 
main non-zero, e.g., contributing to the rms velocity in 
this case (see H2.7D . The product of fluctuations in veloc¬ 
ity and temperature give rise to the enthalpy flux; those 
in velocity and composition give rise to the composition 
flux. 

A stellar evolution code must step over the shorter 
turnover time scales (weather) to solve for the evolution¬ 
ary times (climate). How can this be done? It requires an 
average over active and inactive cells. The steady-state 
limit of the Lorenz equation seems to give a reasonable 
approximation to its average behavior, filterin g out the 
chaotic fluctuations (1 Arnett fe Meakinll2011bD . Instead 
of du/dt = 0, we use 


du/dt = du/dt + (u • V)u, 
-A (u • V)u. 


(17) 


We apply the same approximation (Eq. [T7]) to Eq. [7] for 
slow stages of stellar evolution. This allows non-local 
behavior, will prove important for our discussion of con¬ 
vective boundaries later in (J3J and can represent ram 
pressure (Reynolds stress) and the flux of turb ulent ki¬ 
netic energy; see also §3.2 in lPorter fc Woodwardl (l2000t) . 
for a discussion of ram pressure in 3D simulations relative 
to MLT. 

Now we have established connections between an ac¬ 
celeration equation (Eq. 0 and (1) MLT, (2) histori¬ 
cal at tempts to exte nd MLT, (3) modern research on 
RTI (lAb arz bH l201f)l). _{4) the impor t ant a dvances of 
IKolmogorovI ( 1941L I1962D and iLorerO (|1963f) . and (5) a 
rational way to step over fluctuations for stellar evolu¬ 
tion. 


2.7. Turbulent Kinetic Energy Equation 

A more rigorous alternative is to use the Reynolds- 
averaged Navier-Stokes (RANS) approach, which di¬ 
rectly averages the fluctuation s over space and time. This 
has been explored by Ca n uto d(.’anut oll20 11 al HEIdl j). se e 
also iXiong. et al.1 (|1997f ) ; iDeng. Xiong, fc Chanl ) 20061) ; 
a detailed comparison with their work, while desirable, 
is beyond the scope of this paper. Canuto uses simula¬ 
tions and experiments from geophysics to effect a closure 
of the RANS equations, while in contrast, our closure of 
the RANS is based on our 3D simulations. 

The turbulent kinetic energy equation (TKE) is ob¬ 
tained by a Reynolds decomposition of the velocity, den¬ 
sit y, and pressure (detailed discussion may be found 
in iMeakin fe Arnett! I2007bl: lArnett. Meakin. &; Yound 

120091 : IViallet. et al.l 120131 iMocak. et al.l 120141 ) In orin- 
ciple the TKE is exact; errors arise from closure, 
i.e., our analytical approximations to the terms in the 
RANS equations are at fault. Well-resolved 3D ILES 
simu l ations show e xcellent agreement with the TKE 
(jViallet. et alJl2013l) . and allow the dominant terms to be 
identified. Being more general than the simpler approx- 
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imations discussed above, the TKE allows us to identify 
and quantify neglected terms. Most importantly, it al¬ 
lows an enormous simplification and compaction of the 
3D numerical data, while that data in turn allows a clo¬ 
sure of the RANS procedure. 

The TKE may be written as dMeakin fc Arnettl 
I2007bf) : 


What approximations would be necessary to make the 
TKE equation equivalent to MLT? In MLT, (1) the net 
flux of turbulent kinetic energy Fk is defined to be zero 
by symmetry, (2) pressure fluctuations are ignored so 
the acoustic flux Fp and pressure dilatation (P'V • u') 
are zero, and (3) the damping length i d is taken to be an 
arbitrary adjustable parameter. Enforcing these gives 


d t {pE K ) + S7-{pE K u 0 ) = 

-V-(F p + F k ) + (P'V ■ u') 

+ {p' g • u') - pe dl (18) 

We use ( q ) and q to denote angular and time averages of 
a quantity q. Primes refer to fluctuating quantities; for 
example u = uo + u', and (u) = uo, and similarly for 
the time average. The turbulent kinetic energy per unit 
mass is Ek = |(u' • u'), a measure of the rms turbulent 
velocity. The acoustic and turbulent kinetic fluxes are 
Fp = P'u' and Fk = pEx u'. The dissipation may be 
written as 

£d = (u' • \i’\u'\)/i, (19) 

a form which we identif y with Eq. [TJ the expression of 
I Koln locorovl (1 191 ill 1 962fl : notice that it involves averages 
of powers of the velocity fluctuation, not the instanta¬ 
neous values. 

Using the RANS approach is equivalent to using the 
bottom panel in Fig. [2] rather than the top; it removes 
the fluctuating activity which cancels (has no net effect), 
while keeping what does not cancel. 

To better understand the implications of the TKE, con¬ 
sider (1) a steady state ( d t (pEK) = 0) with (2) no back¬ 
ground motion (uo = 0). Then the TKE reduces to the 
divergence of the fluxes V-(Fp + Fk), balancing the net 
result of two source terms (P'V • u') and (p'g • u'), and 
a damping term —pe d : 

V-(F P +F k ) = (P'V ■ u') 

+ {p' g • u'} - pe d . (20) 

This may be integrated over the convection zone (taking 
the surface fluxes to be zero or small at the boundaries), 
and if we ignore the pressure dilatation (P'V • u' for the 
moment, gives an expression for the damping length i d . 


td = 



(—g • u' 

cz Po 


dm. 


( 21 ) 


which is a global condition that must be satisfied to 
be consistent with Kolmogorov damping, which also re¬ 
quires that id is approximately the depth of the turbulent 
region. This characteristic length scale is a fundamental 
property of turbulence, and is generated robustly in the 
numerical simulations. 

Eq. [Til mi ght be rega rded as a generalization of the 
iRoxburghl (IlfiM Il992h integral constraint to include 
damping by turbulence. Notice that i d . which appears in 
both Eq. [3 and Eq. [2U must be solved for consistently; 
it tends to be a slowly-varying function, of order of the 
convective zone depth. Eq. [2T] involves some of tire im¬ 
portant “bulk” properties discussed bv iCanutol (11992D . 
and is a statement of a global balance between driving 
and damping. 


(p'g ■ u') = pe d . (22) 

This is the local version of the global balance in Eq. I2T1 
it is equivalent to the Bohm-Vitense cubic equation 
of MLT for the appr opriate choice of mixing length 
(jSmith fc Arnettl 120141) . 

This approximation leads to a series of errors: (1) 
Symmetry between up-flows and down-flows is bro¬ 
ken by stratification, so that turbulent kinetic energy 
fluxes are not general l y zero (iStein fe Nordlundl 119891 
iCattaneo. et al. 1 119911 [Canutolfl992D . This is a quali¬ 
tative error. (2) Pressure fluctuations may not be ig¬ 
nored for strongly stratified convection zones. This is a 
quantitative error. IViallet. et all ( 2013 )1 find that accel¬ 
eration by the pressure dilatation term is comparable to 
that from buoyancy. (3) The d amping length may not be 
freely adjusted if the relation of lKolmogorovl (11 91 ill 1962l ) 
is to be satisfied. Such adjustments are usually necessary 
to compensate for a lack of non-locality in atmospheres 
due to the lack of ram pressure, and deeper into interiors 
due to a lack of kinetic energy flux (the two parameters 
discussed in regard to 3D atmospheres in 1 12.31) . 

2.8. The lPasetto. et~al\ H201A ) model 

Our efforts have been three-fold: (1) construction of 
accurate numerical solutions of the Navier-Stokes equa¬ 
tions which exhibit turbulence, (2) theoretical analysis 
of these solutions in the RANS framework to determine 
the most important features, and (3) invention of simpler 
analytic representations which cap ture the essent i al fea- 
tures of the numerical solutions. iPasetto. et al.1 (|2014f ) 
have presented a novel analytical theory of convection in 
stars which does not contain a mixing-length parameter; 
this is an alternative to (3) above, and it is of interest 
to compare how well it agrees with both our numerical 
solutions (1 and 2), and o ur an alytic approximations (3). 

As we have shown in m the natural length scale 
for convection is the dissipation length for the turbu- 
lent cascade. Part of the foundation of the model of 
IPasetto. et all (|2014 ) is the use of potential flo w and the 
Bernoulli equation (lLandau fe Lifshitzl JI959 1 ), Eq. 10.7 
in §10), which result from the Euler equation, not the 
Navier-Stokes equation. Their theory seems to be equiv¬ 
alent to assuming the process occurs on a scale much less 
than the size of the convective region, so that there is no 
way to define a length scale for turbulent dissipation. In 
contrast, following Kolmogorov 1 42.II) . the length scale 
in our theory is the size of the turbulent region, which 
is not arbitrary but determined by the turbulent flow. 
Our length scale is not an assumption (as in MLT) but a 
consistent and robust result of our simulations. It is the 
length scale over which driving and damping of turbu¬ 
lence balance ( 42.31) . In order to describe the turbulent 
cascade, a complete theory must deal with the whole tur¬ 
bulent region. 

Is the theory of IPasetto. et al.1 (|2014T ) physically cor- 
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rect? Stellar interior convection is extremely turbulent, 
so the question becomes: w hat are the errors intro¬ 
duced by ignoring turbulence? lLandau fe Lifshitzl (|1959f l 
give a careful discussion of the applicability of poten¬ 
tial flow (their §9), and they note that the validity of 
Bernoulli’s equation is limited because of the formation 
of boundary lay ers i n which viscous effect s must be in¬ 
cluded (see also iPrandtl fe Tietiensl 119341) . Stars have 
large Rey nolds numbers, so that tur bulent boundary lay¬ 
ers form (jLandau fc Lifshit^ (1 9591 . Chap. Ill), as they 
do in our simulations (Fig. [3]). The Pasetto theory, like 
MLT, ignores boundary layers and turbulence, as well as 
composition gradients. 

A basic assumption of the iPasetto. et al.l (120141) the¬ 
ory is that velocities of lateral expansion are much larger 
than those of the vertical rise of convective elements 
(their §4.2). However, the simulations show average ve¬ 
locities in the turbulent region which are not strongly 
biased t oward the laterial directio ns: this was already 
clear in iMeakin fe Arnettl (l2007bl) . (their Fig. 6), and 
has held t rue for subsequent simulations with refine d 
resolution (|Viallet. et all 120131 ICampbell. et al~l 120151) . 
The rms velocity in the radial direction is actually 
(arqerthan the lateral rms veloci ties, rather than smaller 
(|Arnett. Meakin. fc Youndl2009D. _ 

A key test presented in IPasetto. et al.l (120141) of their 
theory is a comparison with MLT0 at r = O.98f?0, well 
inside the super-adiabatic region (SAR) at r ~ 0.9985i?g 
in the Sun. It is the inefficient convection in the SAR 
which determines the solar radius in calibrations of stel¬ 
lar evolutionar y codes, so that a test in the SAR would 
be instructive. IPasetto. et al.l (12014 ) state “Convective 
elements in this region have low thermal capacity, so that 
the super-adiabatic approximation can no longer be ap¬ 
plied, and the temperature gradient of the elements and 
surrounding medium must be determined separately”. 
The theory in its present form may not yet be applicable 
to the SAR. 

The value of the Pasetto theory may prove to lie in its 
significant conceptual differences from MLT, and in its 
use as a null case to provide insight into the effects of 
turbulence. 

3. BOUNDARIES AND BOUNDARY LAYERS 

It has been assumed that because deep convection is 
adiabatic, MLT may be used without problem for stan¬ 
dard stellar evolution in deep interiors. This ignores the 
effects of the velocity field. Realistic boundary physics 
requires more than the adiabatic assumption; it requires 
dynamics to define the b oundary, and hence the size 
of the convective regions (lArnettll 19911 lAsida fc Arnettl 
l2000t IMeakin fc Ariiettll2f)()6l I2007H) . 

Because, unlike MLT, Eq. [7] and its variants have a 
spatial derivative, the edges of the convective zones may 
he found by simply integrating the acceleration equation 
to find the zeros of the velocity. 

In this section we begin by discussing several issues re¬ 
lated to boundaries. We stress the importance of Peclet 
number variation (» We critically review current 
practice regarding artificial diffusion, real diffusion, semi¬ 
convection, and imposed boundary criteria (icinrao 
13.51) . Then we discuss the similarities and differences be- 

18 As our title suggests, we attempt to go beyond MLT. 


tween convection in stellar atmospheres and deep inte¬ 
riors 1 43.61) . In 43.71 we present new numerical results 
concerning convective boundaries (the development of 
braking regions, which do not appear in MLT). In 43.81 
we then analyze these results, showing that they emerge 
from simple considerations of physics, which may be used 
to construct approximations for use in stellar evolution¬ 
ary codes. 

3.1. Peclet number: radiative diffusion 

For the oxygen-burning shell, the temperature T has 
an abrupt jump inside the mixing region (radius r ~ 
4.3 x 10 8 cm in Fig. H). Pressure is continuous through 
the boundary containing this transition, so that the 
density curve ha s a corr esponding dip; see Fig. 2 in 
IMeakin fc Arnettl (|2007bD or Fig. 5 in iViallet, et all 
( 20131) . This implies a steep increase in entropy; as evo¬ 
lution continues this entropy jump grows, and the tran¬ 
sition region narrows. Such steep gradients in T are a 
consequence of cooling by neutrinos. They are not seen 
in earlier, photon-cooled stages of evolution and can only 
be supported for times short compared with timescales 
for thermal diffusion and electron heat conduction. This 
is easily the case for oxygen burning because of high opac¬ 
ity and short evolutionary times (~ 10 5 sec). 

The Peclet number is defined as the ratio of the advec- 
tive transport rate to the diffusive transport rate of the 
physical quantity being transported, which here we take 
to be thermal energy, so 

thermal advection rate 

Pe = -;-;-. 

thermal diffusion rate 

In oxygen burning, radiative diffusion is slow while ad¬ 
vection occurs rapidly, giving large Peclet numbers (for¬ 
mally infinite since radiative diffusion was small enough 
to be neglected in some simulations; the infinity results 
from the denominator in the definition being a negligible 
term, not from any exceptional behavior of the physics). 

This contrasts with the situation in stellar atmo¬ 
spheres, in which the radiative diffusion becomes faster 
than advective transport, so that Pe < 1. This dif¬ 
ference in Peclet numbers suggests the possibility of a 
fundamental flaw in the notion that observations of stel¬ 
lar atmospheres may be sufficient to define the n a ture o f 
deep stellar convect ion. See discussion in iZahnl (119911) ; 
lYiallet. et al.l (|2015l ). 

3.2. Artificial diffusion 

Peter Eggleton took an early step in dealing with steep 
gradients in composition, with the introductio n of a d if¬ 
fusion operator which he stressed was ad-hoc (lEggletonl 
I1973D . This numerically advantageous procedure has 
been widely adopted for stellar evolution, even though 
it has the potentially worrisome mathematical property 
that it increases the order of t he spati al derivatives in the 
equations to be solved. The lEgglctoiJ (I1973D equation is 



where X is the mass fraction, m is the lagrangian mass 
coordinate, a = vml^ml (47rr 2 p ) 2 is the effective diffu¬ 
sion coe fficient, and 1Z is the nuclear reaction network 
matrix (lArnettl I1996D . This is equivalent to modeling 
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convection as “turbulent diffusion.” The left-hand side 
is the heuristic diffusion operator; the right hand side 
is the reaction network operator. The actual composi¬ 
tion flux is related to the co-mov ing derivative on the 
right-hand side; see lArnettl 119961 §4.6. Eggleton inte¬ 
grates over the convection zone to eliminate that spatial 
derivative; usually it is simply ignored in stellar codes. 

The Eggleton approach is equivalent to approximating 
the composition flux 

F y = pA(u'Y') (24) 

by a “down-g radient” expression (critically discussed by 
lCanutolll992l'l . 

Fy —> pA u{—IdY/dr). (25) 

Direct comparison with simulations shows that this can 
be qualitatively wrong (b y two orders of magnitude ). For 
a contact discontinuity (ILandau fc Lifshita I1959L §81), 
Fy —> pAuAY, as in Eq.[24j not pAu(—I/Ar)AY —)■ oo, 
as in Eq. [25] Proper scaling requires that I —>• A r at a 
boundary if Eq. [25] is used. 

As Eggleton intended, the algorithm smooths steep 
gradients, but sometimes faster than real physical pro¬ 
cesses do, as Eggleton warned. To the extent that gra¬ 
dients in abundance need to be correctly represented 
(e.g., for ionic diffusion, or density structure), the down- 
gradient approximation (in Eq. [33] and Eq. 1351) . is 
questionable. I n particular, fluxes dire c tly computed 
in si mulations dMeakin fc Arnettl l2007bt lYiallet. et al.l 
l2013f) show that th e down-gradient ap proximation fails 
in boundary layers dMocak, et alj|2014h . 


3.3. Ionic diffusion 

While real atomic (ionic) diffusion is thought to be slow 
in stars, the diffusion operator is second order in space 
derivatives, so that it becomes important in steep com¬ 
position gradients, i.e., boundaries. Georges Michaud 
has led in the application of true diffusio n processes and 
radiative levitation to stellar evolution (IMichaudlll970L 
119911 : [Michaud. Richer, fc Richardll2007l) . Recently these 
proce sses have been applied to horizontal branch and sdB 
stars (IMichaud. Richer, fc Richardll2005l l2007t IHu et al. ' 


2008 1 120091 1201(1 1201 It IMichaud. Richer, fc Richard 


2011 Bloemen et al.1 I2014h . Gravit ational settling 
( Hu et al.H2009f) andradiative levitation (IHu et al.ll 201 ll 
are important to ( 1 ) recover t he iron - group opacity bu mp 
that excites the pulsations (ICharninet et al.1 11997T) in 
those stars, ( 2 ) obtain the correct position of the in¬ 
stability strip in the log <7 — T e ff diagram, and (3) help 
in u nderstandin g their observed atm ospheric abundances 
([Michaud. Richer, fe Richardll2011lh 
Because the lEggletonl ( 1973T) diffusion uses a difference 
operator similar to that for ionic diffusion (second order 
in space), and may reduce the gradients which drive that 
diffusion, care should be taken that the algorithmic dif¬ 
fusion does not cause errors in the real diffusion (e.g., see 
ISchindler, Green fc Arnettl 120151) . 


3.4. Semi-convection 

In stellar physics, the idea of s emi-convection 
has s pawned various algorithms (e.g.. iSchwarzschild 
19581 IStothersI 119631: ICastellani. Giannone. fc Renzini 

ToTiaHhT iDemaraue fc Mengell Il972t iSweigart fc Grossl 


119761 : iDorman fc Roodll 19931) . some of which seem to be 
physically and numerically inconsistent with others. The 
term “semi-convection” refers to a mixing process which 
occurs in a region that is stable according to the Ledoux 
criterion but unstable according to the Schwarzschild cri¬ 
terion. It generally is thought to involve mixing of com¬ 
position, but not significant enthalpy. The composition 
profile may be adjusted to marginal stability according 
to the Ledoux criterion. 

Semi-convection is also often discussed as a dou¬ 
ble diffusive instability, involving an int eraction be¬ 
tween radiative diffusio n and ionic diffusion (|Snrui mm 
lLattanzio. et al.112014T) . Although both radiative and 
ionic diffusion may be included in a ID stellar code, 
this does not capture their interaction and 3D dynamics. 
Semi-convection may be related to oceanic phenomena 
(thermohaline mixing) in which heat flow and salt con¬ 
centration play the doubly-diffusive roles, a nd w hic h have 
a long his t ory of study (e.g., see Ch a p. 8 inlTurneillT973l : 
lGilMl982lh iRosenblum. et al.l (|2011lkfWood. et al.11)201311 
give an extensive discussion with numerical simulations 
based on the oceanic model, and conclude that, while the 
problem can be solved in the planetary range of parame¬ 
ter space, the stellar case requires a large extrapolation. 
This difficulty may be further exacerbated by the indi¬ 
cation that many such regions in stars are bath ed in a 
flux of g-mode waves dMeakin fc Arnettl l2007bf) . which 
are a nonlocal eff ect that may comp licate the analysis in 
a nontrivial way dMocak. et al.lf2014h . 

Even with these uncertainties, there are energetic con¬ 
straints (see Eq. [301) which must be obeyed. The amount 
of mixing possible is limited by the energy available to 
mix, which is generally taken to be related to the excess 
V r — V 0 , so that luminosity is used to supply the energy 
required to mix. 

3.5. Imposed Boundaries 

MLT, as a local theory, must be supplemented by ad¬ 
ditional assumptions ab out behavior at the boundaries 
of the convection zone dSpiegell I1971LI1972T ) . These are 
usually discussed in terms of linear stability theory, i.e., 
in t erms of the Ledoux and the Schwarzschild crite¬ 
ria dKippenhahn fc Weigertl fT990f) being positive. The 
Schwarzschild criterion for convective instability is de¬ 
fined by 

<S = V r - V 0 . (26) 

Here V r is what the dimensionless temperature gradi¬ 
ent would be if all the luminosity were carried by radia¬ 
tive diffusion and V a is the adiabatic gradient (see Ap¬ 
pendix). The Ledoux criterion for convective instability 
has a composition dependence, and is defined by 

By 

C = V r - V a - j-Vy. (27) 

PT 

The last term is written as by 

iKippenhahn & Weigertl (jl990l h §6.1, their Eq. 6.12. 
The B factors are defined as in H2.4I above. Notice that 
positive 4^-Vy and positive both inhibit mixing. 

Neither of these choices seems satisfactory. They have 
no dependence upon the vigor of the flow on the unstable 
side of the boundary, which clearly must make a differ¬ 
ence. 
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Linear perturbation theory examines the instability 
of a stable region, treating both sides of the bound¬ 
ary equally. In reality they differ: one side is convec¬ 
tive. The stiffness of the non-convective side is mea¬ 
sured by th e Brunt-Vaisala (buoyancy) fr equency N, (see 
Eq. 6.18 inlKipoenhahn fc Weigertll 199(1 and Eq. 3.73 in 
lAerts, et alJl201Clh. where 

^ 2 = -^( V e- V a-(/W/?T)W) = -^£. (28) 

N is the frequency of elastic rebound from a perturba¬ 
tion; it is imaginary in convective regions. Here V e is 
the dimensionless temperature gradient relevant^ to the 
perturbed element. On the non-convective side of the 
boundary, it may be the same as V r above, giving the 
second equality, which refers to the tendency to restore 
stability in the radiative region. 

A delicate point i s the value of Vy near the bound¬ 
ary (iGabriel. et al.l I2014D . By what mechanism does 
mixing occur? What is the structure of the partially 
mixed region of transition between well-mixed and un¬ 
mixed? Present practice in stellar evolution is to use the 
Schwarzschild criterion, which has no Vy, so that these 
issues may be ignored, or to use the Ledoux criterion 
with one of the prescriptions for semi-convective mixing 
(see H3.4I) . 

Such interfacial issues have long been studied in the 
fluid dynamics and geophysics communities; see iTurnerl 
dl973l) for an extensive discussion. The Richardson num¬ 
ber is defined as some measure of 

potential energy needed to mix 
1 ~ kinetic energy available to mix ' 

The linear condition for ability of a layer to resist shear 
is the “gradient” Richardson number Ri. 

Ri = N 2 /(du/dr) 2 > 1 (29) 

is stable; larger stiffness ( N 2 ) and less swirling 
(( du/dr ) 2 ) ten d toward stability . In t heir discussion 
of entrainment. iMeakin fe Arnettl (l2007bl) used a “bulk” 
(i.e., non-local and non-linear) Richardson number which 
involved an integral over the region around the boundary. 

In the absence of g lobal rotation, a layer having con¬ 
stant total entropjo is energetically neutral with regard 
to mixing. If after a mixing episode, the luminosity re¬ 
turns to its value for radiative balance (V r is unchanged), 
then the additional energjEH required to remove the sta¬ 
ble compositional stratification is 

E mix = gH P p T {C ~ S) = gH P {-/3 Y V Y ). (30) 

Both /3t and j3y are intrinsically negative in stars. 
If this energy E m i X changes sign, mixing may oc¬ 
cur which is driven by the gradient in composition 

19 The exact meaning depends upon the assum ed flow , and i s dif¬ 
ferent for MLT a nd th e Lorenz model fsee lArnett fc Meakinl2011bl ; 
Smith & Arnett 2014; and Table[l}. 

20 See, e.g.. I Arnettl (1 19961 1 for explicit derivations of all compo¬ 
nents of the entropy (Appendix B), and of the total energy of the 
star (Appendix C). 

21 This is the change in internal energy due to composition 
change, keeping temperature and pressure constant. 


(jMocak. Meakin, fc Miillerll2011bf) . Using a specific ki¬ 
netic energy of \u 2 : a Richardson number may be con¬ 
structed, 

Ri = 2gH P {-l3 Y V Y )/u 2 . (31) 

Here the traditional Ri > 1/4 is a plausible condition for 
stability, at least roughly. 

3.6. Solar convection 

In their pioneering work on solar convection, 
IStein fe Nordlundl (119891) carefully explored the topology 
of convective flow below the photosphere: converging, 
cool downdrafts being dominant, with radiative cool¬ 
ing provi ding the entropy deficit w hic h driv es the cir¬ 
culation. iFrevtag. Ludwig, fe Steffen] (I1996T) examined 
shallow (weakly stratified) convection, driven by atmo¬ 
spheric cooling, and emphasized the importance of the 
atmosphere in determining the nature o f the convec¬ 
tion zone. As deep interior convection (|Arnettl 11994 
iBazan fe Arnettl 119941 has no atmosphere, atmospheric 
physics can have no strong role there (the circulation is 
driven by nuclear burning). Furthermore, the bottom 
boundary, which could b e ignored in the simulations of 
IStein fe Nordlundl (I1989T ) . may be important for the de¬ 
tai led effects of solar c onvection on the interior. 

iSchwarzschildl (119581) . §11, showed that, for stellar inte¬ 
rior models, the atmosphere could be represented by an 
entropy jump between the photosphere and the adiabatic 
(deep) convective region. This entropy jump is a primary 
parameter for determining the depth of the convection 
zone. The atmospheric model is crucial for predicting 
spectral features for a given entropy jum p, but has a 
weak influence on that entropy jump itself ( Tanner, et al 

l20lll2Ml . 

Many features of the atmospheric and deep inte¬ 
rior simulations are similar, leading to the idea that 
atmosph eric physics, however cru c ial for spectral for- 
mation (|Stein fc Nordlundl 119981 : iMagic. et al.l I2013L 
HOUD, may be treated as a boundary condition issue 
rathe r than a key fea t ure o f deep turbulent convec¬ 
tion. IMeakin fe Arnettl (1201(1 ) showed that the general 
characteristics of the flow in solar convection (narrow, 
fast down-flows with broad, s low up-flows and accel¬ 
eration by pressure dilatation, IStein fe Nordlundl 119891 : 
IViallct. et al.l 120131) . require only localized top cool¬ 
ing and stratification. Global simulations of the so¬ 
lar convection zone are necessarily less well resolved for 
comparable computa tional resources; the simulations of 
iMiesch. et al.l (I2007D are beginning to show turbulence, 
but may require finer zonin g to deal wi t h some detail s of 


the turbulent flow (e.g., 

Hanasoge, Duval & Sreenivasan 

2012: Brandenburg 2015 

)■ 


3.7. Deep interior convection 

The simplest of stellar convection zones are cooled by 
the local processes (cooling by neutrino emission and 
heating by nuclear burning), rather than the non-local 
processes (radiative transfer), giving a cleaner example 
of the dynamics of boundaries for deep convection. A 
slightly more complex case is a conv ection zone wi th heat 
conduction by radiative diffusion; IViallet. et alJ (|2013l ) 
consider both. These two cases cover almost all of the 
conditions relevant to stellar evolution, except the outer 
layers simulated in 3D atmospheres. 
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Table 2 

Integral Properties of Convection Zone Regions 


variable 

symbol 

total (CZ+BL+BL) 

lower BL 

upper BL 

mass 

Attt/Mq 

0.9205 

0.0161 

0.1150 

depth 

Ar/10 s cm 

4.460 

0.078 

0.587 

kinetic energy 

KE/10 46 erg 

8.608 

0.255 

0.561 

buoyancy luminosity 

Lbuoy/10 45 erg/s 

4.576 

-0.0342 

-0.0492 

pressure 

AlnP 

2.032 

0.046 

0.228 

number of zones 0 

A i 

236 

8 

23 


a The total number of zones in the radi al direction was 400 in Meakin & Arnett 2007b (this table, medium 
resolution), 800 in IViallet. et al.l 120131 (high resolution), and 1536 in Fig. [T] The basic features appear 
even at lower resolutions. 


For the oxygen-burning shell, some integral properties 
of the main convective region and the braking layers are 
summarized in Table [2j About 14 percent of the mass 
and 15 percent of the thickness of the total convection 
zone are in the boundary layers (upper BL and lower BL), 
as is 8.5 percent of the turbulent kinetic energy. These 
boundary regions provide deceleration (braking) of the 
vertically directed flow, allowing it to remain bounded 
by the convective volume. If the buoyancy flux is q = 
—g{u' z p')/po, then the rate at which turbulent kinetic 
energy increases due to buoyancy in a region a, is 



which is positive in the middle region, but negative in the 
boundary regions. These regions of negative buoyancy 
are a robust qualitative feature of the s imul ations, datin g 
back to early 2D work dHurlburt. Toomre. fe Massaguerl 
11984 lArnettl I1994D . In the oxygen-burning shell they 
reduce the driving of turbulent kinetic energy by only 
1.8 percent. 

Table [2] shows the depth of each region in pressure 
scale heights (A Ini 3 ). The depth of the boundary zones 
is not a universal constant in In P, but varies by a factor 
of 5 between top and bottom. The last line gives the 
number of zones in each r egion for “medium” resolution 
dMeakin fe Arnettll2007bD ; the lower boundary region is 
most demanding, having a steep transition from convec¬ 
tive to stable stratification. 

Little of the kinetic energy is lost in the boundary re¬ 
gions, so Lbuoy provides a good first estimate of the rate 
of generation of turbulent kinetic energy. These regions 
contain 17% of the mass in the “convection zone”; most 
of this comes from the upper layer, which has less ex¬ 
treme stratification. 

Fig. [3] shows the buoyancy flux versus radius, averaged 
over 100 seconds, for the oxygen-burni ng she ll sim u lation 
(OB); more detail may be found i n dMeakin fc Arnett 
200m; I Arnett. Meakin. fc Younel 120091 IViallet. et al. 
20131) . The buoyancy flux —u • gp'/po, is the rate of 
work done by gravity (IZahnlll99lD . It is the rate of flow 
of buoyancy, —gp'/po, and has units of energy per unit 
mass per unit time (e.g., erg/g/s). Over most of the 
convective region it is proportio nal to the enthalpy flux 
IjArnett. Meakin. fc Youngj[2009 ). 

Fig.[3]shows that the convective zone simulation is nat¬ 
urally split into three regions, separated by two bound¬ 
aries. The regions above and below are stable. The mid¬ 
dle region is relatively uninfluenced by the boundaries; 



Figure 3. Buoyancy Braking averaged over 100 seconds (~ 
2 transi t times) at shell boundar i es for oxygen burn ing: q versus 
radius (Meakin & Arnett 2 007H ; IViallet. et all 12013 ) . The buoy¬ 
ant acceleration changes sign near the boundaries of the convection 
zone, giving braking rather than positive acceleration. 

it is characterized by positive fluxes of buoyancy and of 
enthalpy, that is, a positive “superadiabatic gradient” 
AV. It is convectively unstable according to both the 
Schwarzschild and the Ledoux criteria. With an appro¬ 
priately choice of mixing length, this middle region can 
be reasonably well approximated by MLT. 

MLT works poorly for the bottom and top boundary 
layers, which have negative values of AV. While the cen¬ 
tral region is defined by positive buoyancy, and positive 
enthalpy flux, outside the convective zone these quanti¬ 
ties are zero, and in the boundaries they are negative. 
In MLT this is impossible because it would imply that 
the velocity in Eq. [5] is imaginary, but in Eq. [7] merely 
implies buoyancy braking, hence the labels “braking” in 
F ig. [3] 

IZahnl (11991D has summarizetS the issue of negative 
buoyancy and con vective flu x in connection with pene¬ 
trative convection. ISchmitt. Rosncr fc Bohn ( 19841 ) have 
discussed the overshoot at the bottom of the solar con¬ 
vection zone in the co ntext of convective pl umes and 
magnetic dynamos, and lSniegel fc Zahnl (119921) have dis¬ 
cussed this in the context of solar rotation and the 
tachocline. In stellar evolution theory (i.e., MLT) the 
existence of these braking regions is obscured by use 

22 See £|2.7l and Eq. I21l for an explanation of “appropriate.” 

23 Compare his Fig. 1 to the right braking layer in our Fig. [3] 
this is a nice prediction of some of the features later revealed in 3D 
simulations. 
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of the Schwarzschild (or Ledoux) linear stability cri¬ 
terion. These braking layers are relate d to issues of 
overshoot and penetrative convection (IVeronisI H963I: 
Massaeruer et al.1 11984 iHurlburt. Toornrc. fc Massaguerl 
198 (tT ) . The braking layers are not a part of MLT but, as 
we shall see (1 13.81) . arise naturally from Eq. [7] 

Fig. H] shows the inner braking zone (the region of neg¬ 
ative buoyancy work) at r ~ (0.433 to 0.445 x 10 9 cm). 
The “hi-rcs” case of lYiallet. et al.ll2013l (768 x 512 2 zones) 
and a still higher-resolution case of iCampbell. et al.l 
120151 (1536 x 1024 2 zones) are shown. In comparison with 
Fig. m the negative “spike” is now well-resolved. A de¬ 
tailed analysis of these simulations will appear elsewhere. 
The degree of numerical convergence is promising, and 
we conclude that such braking zones are a robust fea¬ 
ture of well-resolved simulations of neutrino-cooled stel¬ 
lar convection. 



Radius (10 9 cm) 

Figure 4. Time-averaged buoyancy work (weighted by a factor of 
4irr 2 p) at lower shell boundary for oxygen burning, versus radius. 
This shows the “hi-res” case of ( IViallet. et alJl2013f > (768 x 512 2 ) 
and a higher-resolution case (1536 X 1024 2 ). The braking zone 
is indicated by negative buoyancy work at (0.433 to 0.445 X 10 9 
cm). Compare to Fig. [3] which shows both the upper and lower 
boundary for the ” medium-res” case. There is a steady convergence 
toward a common asymptote as resolution increases, and the two 
cases shown here are virtually identical, except for small variations 
in averaging due to differences in time step size. 


The radial velocity becomes small in the braking re¬ 
gion, while the transverse velocity extends deeper be¬ 
fore it also becomes small. The convective motion turns, 
and a small (mostly g-mode) wave velocity remains. The 
composition gradient is steeper than would be predicted 
by algorithmic diffusion (Eq. l23ll , and begins at the bot¬ 
tom of the braking region. The boundary composition 
profiles are smooth and self-similar when time-averaged. 
This suggests that the turbulent spectrum has a consis¬ 
tent net effect on the composition profiles and on the 
mixing, and therefore this interface should be amenable 
to approximation over time-steps in ID evolutionary cal¬ 
culations. 

For oxygen burning, the composition gradient in the 
boundary layer is not well-represented by conventional 
turbulent diffusion theory which requires a span of many 
“turbulence mean-free-paths” per density scale height 


(jAmsden fc Harlow! I1968D for validitjH In MLT, the 
span is a fraction of a scale height (see A In P in Table [2j 
for oxygen burning. The small length scales are accom¬ 
panied by small time scales for change, so that a steady 
state model may be appropriate. 

3.8. Dynamics and Braking Layers 

Fluid motion in a star m ay be separated into tw o 
fundamentally different flows dLandau fc Lifshit3ll959 )) : 
solenoidal flow (divergence free: V-pu = 0) and poten¬ 
tial flow (curl free: V x pu = 0), which together repre¬ 
sent the Helmholtz decomposition of an arbitrary vector 
field. Potential flow is associated with wave motion and 
solenoidal flow (vorticity) is a feature of turbulence. A 
striking separation in the nature of the flow is visible 
at boundaries between these types of flow; see the dis ¬ 
cussion of boundar y layer s in lPrandtl fc Tietiensl (119341): 
lLandau fc Lifshit3 (119591 ) . and Fig. 19 in IViallet. et al l 
( 20131) . This separation in types of flow is c losely 
related to wave generation and propagation ( Press| 
19811: iPress fc Rvbickil 119811: Goldreich fc Kurnarl 19901 : 

Goldreich. Murray, fc Kuniarl 1994D. 

The structure and nature of these boundary lay¬ 
ers is important for estimation of the rate at which 
turbulent flow moves into or from non-turbulent 
regions—the growth and recession of convective zones. 
iMeakin fc Arnett : (]200711 ) had about 8 zones across the 
lower boundary laye r for “medium” r esolut ion; see also 
iHerwig. et al.l (|2014D . IViallet. et al.1 (I2013D had double 
the resolution across the convective zone (twice as many 
radial points), but the boundary layer became physically 
narrower. Rece nt simulations at still h igher resolution 
(see Fig. 0] and ICampbell. et al. I I2015D show that the 
lower boundary layer has about 20 zones and the same 
physical depth. The computed entrainment rate may be 
affected by numerical viscosity, so that lower resolution 
simulations will give overesti mates. _ 

The “medium” resolution of lMcakin fc Arnettl (l2007bD 
was sufficient to give numerical viscosity (Reynolds num¬ 
ber) similar to that of laboratory experiments on entrain¬ 
ment, but not of stars. Coarse resolution in those sim¬ 
ulations may have been a par tial cause of the difficul¬ 
ties found by iStaritsinl (120131) in an attempt to a pply 
the entrainment rates of IMeakin fc Arnettl (|2007bl) for 
oxygen burning directly to main sequence stars. The 
real entrainment rates for stars should be smaller. An¬ 
other issue is that oxygen burning and h ydrogen b urning 
have very different Peclet numbers (IViallet, et al .1120151) , 
which can affect the entrainment rate (see below). 

Here we construct a simple but dynamically consistent 
picture of a convective boundary. This is illustrated in 
Fig.0 which shows the driving, turning, shear and stable 
regions. At its most elemental level, the velocity vector 
must turn at boundaries; that is, flow must turn back to 
stay inside the convective region. We do not assume that 
“blobs” disappear (like MLT). Most of the momentum is 
contained in the largest scales, so we focus on the average 
dynamics at these scales, and the simplest flow patterns. 

The magnitude of the acceleration required to turn the 
flow is just the centrifugal value u 2 /b where b is the radius 
of the turning region and u the relevant velocity. Using 

24 The problem is similar to that in a stellar photosphere, in 
which radiative diffusion must give way to radiative transfer. 
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Figure 5. Simplified schematic of a convective boundary. The 
length b corresponds to the radius of curvature needed to reverse 
(contain) the flow ( u r —> —u r ). The centrifugal acceleration is pro¬ 
vided by pressure fluctuations (see text). The boundaries oscillate 
due to surface waves. The radial direction is denoted by r and the 
transverse by h. Orientation is for the top of a convection zone; 
the bottom may be described by appropriate reversals. 

Eq. □ in the steady state limit, and taking b ~ A r -C £, 
the radial component of the acceleration equation be¬ 
comes 

u r du r /dr ~ A(^u 2 )/Ar ~ B, (33) 

where B is the acceleration due to buoyancy and pressure 
fluctuations (Eq. [Gl and HA.21) . So far we have considered 
the top of a convective zone; the bottom of a convection 
zone behaves si milarly if care is taken with signs. _ 

Sim ulations dMeakin fc Arnett] l2007bl : lYiallet. et al.l 
I2013T) show a consistent pattern in velocity and compo¬ 
sition structure in the boundary layers. Moving toward 
the boundary from the interior of the convection zone, we 
find (1) the radial velocity u r decreases, (2) the pressure 
fluctuations P' increase, and (3) the transverse velocity 
Uh increases to a maximum and then decreases, joining 
on to a finite and small rms velocity due to wave motion. 
The transition to small rms velocity occurs at about the 
same point that the composition changes from being well- 
mixed to supporting a radial composition gradient. This 
pattern holds for both top and bottom boundaries. 

The dynamical equations we use are derived in Ap¬ 
pendix [A] We use HA.21 the same quasi-steady state and 
thin shell (b <C £) approximations, and choose an iner¬ 
tial frame in which a hydrostatic background is assumed. 
Near the boundary, the radial component of the acceler¬ 
ation is essentially just 



The buoyancy force (the first term on the RHS) is parallel 
to the gravity vector g, which is radial, and provides 
no transverse acceleration. Baryon conservation implies 
that this reduction in the radial velocity alone will give an 
increase in density (matter accumulates), which gives an 
increase in the pressure fluctuation P' as the boundary is 
approached. The two transverse components of velocity 
satisfy 

Uhduh/dh= --— -dP'/dh. (35) 

Po + P 

The transverse motion requires a transverse accelera¬ 


tion which is provided by a pressure excess (see also 
iStein fe Nordlundl 11989H at the point of contact of the 
plume with the bou n dary (note the similarity to the RTI, 
H2.5I and lSchmitt. Rosner fc Bohnlll984h . 

This same pressure excess also implies a radial accel- 
eration of the boundary , makin g the boundary undulate 
dMeakin fc Arnettll2006l l2007blf . In addition to the hori¬ 
zontal force from the pressure excess, the buoyancy force 
is negative, so the net effect on the flow is to complete 
the turn. The turning region has a width b = r 2 — rq; 
this material is well-mixed because it moves back into the 
convective region after it turns. Thus the region r 2 — rq 
might be termed the “over-shoot” region, and we are dis¬ 
cussing the dynamics of “overshoot”. 

Fig. [4] shows our highest resolution simulation of the 
most demanding boundary; does this simple model of 
boundary dynamics work for it? The orientation is re¬ 
versed for the bottom boundary, so r m i X < r 2 < rq 
in this case. The steep drop in buoyancy work at 
r ~ 0.433 x 10 9 cm corresponds to r m i X and the “shear” 
region in Fig. [5j which can maintain a composition gra¬ 
dient because the velocity is due to wave motion. At the 
radius r 2 , at which the radial component of the velocity 
is u r ~ 0 , the flow is transverse to the radial coordi¬ 
nate (uh ^ 0 ), so there is a shear layer at this surface 
which wi ll be unstable to the K elvin-Hclmholtz (KH) in¬ 
stability (lChandrasekhaiil961f) . The partial mixing layer 
extends to radius r mix (at which Uh ~ 0 ) and contains 
this KH layer. The linear condition for ability of a layer 
to resist shear (stability against mixing) is the “gradi¬ 
ent” Richardson number, Ri > 1/4. The Brunt-Vaisala 
frequency N ~ 3 s -1 is evaluated in the stable region, 
near the boundary, and may be sensitive to resolution. 
The shear velocity is Uh < 0.8 x 10 7 cm/s, and from this 
crude estimate r m i X — r 2 ~ Uh/ZN ~ 10 6 cm. This small 
length is consistent with the steep “cliff” in Fig. [4] 

Both terms in B (Eq. EH) act to turn the flow, and 
are comparable in magnitude. A crude but interesting 
estimate follows if we take B ~ gfdpHp AV, where the 
AV is an average value over r 2 — rq. The turning radius 
in units of local pressure scale height is then 

b/Hp ~ A r/Hp ~ A(iu 2 )/ g(3 T H P AV, (36) 

which is related to the inverse of a Richardson number; 
compare to Eq. [Vjl and [Til Both A(|u 2 ) and AV are 
negative here, giving a positive ratio. The use of Eq. [7] 
automatically leads to an approximate Richardson num¬ 
ber criterion for the edge of the convective region, with¬ 
out the need of an additional imposed boundary condition 
beyond the requirement that u 2 becomes small (see H3.5D . 

The minimum in buoyancy work at r ~ 0.437 x 10 9 cm 
corresponds to rq, the edge of the braking region and the 
“turn” in Fig. [5] At r ~ 0.443 x 10 9 cm the buoyancy 
work becomes positive, so that this corresponds to rq 
and the beginning of the “driving” region, at which AV 
changes sign. Contrary to MLT, the radius rq, at which 
the Schwarzschild criterion is zero, is not at the boundary 
of zero convective motion. 

How does this braking region develop a negative buoy¬ 
ancy? Suppose the region r 2 to rq is well mixed, to 
uniform composition and entropy. There is no braking, 
so convective flow is unabated to the composition gra- 
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dient beginning at r^- Vigorous entrainment erodes the 
boundary, causing a thin layer of partially mixed matter, 
which contains the heavier nuclei from below the oxygen 
burning shell. This makes the buoyancy more negative, 
establishing a braking layer and reducing the rate of en¬ 
trainment. The braking layer grows until the entrain¬ 
ment rate balances the rate of mixing into the edge of the 
convection zone. If the braking layer is too large, such 
mixing will reduce it; there is negative feedback. The 
braking layer is thinner than the convective zone, so the 
time scale is shorter than the turnover time ( H2.1D . and a 
quasi-steady state can be set up. This simplistic analy¬ 
sis (which ignores fluctuations) indicates some of the dy¬ 
namics involved with the braking layers and composition 
boundaries. Furth er analysis with the new higher reso¬ 
lution simulations dCamnbell. et al.~ll2015t iMeakin et al.1 
i2oia is in progress. 

This limiting case (“elastic collision”) is a reasonable 
approximation for the time averaged behav ior of the oxy¬ 
gen burning shell (IMeakin fe Arnettll2007bl) , in which ra¬ 
diative diffusion (and electron heat conduction) are slow; 
here Ttum ~ 0.6 sec, while the radiative diffusion time is 
Tdiff ~ 3 x 10 7 sec. A measure of the heat lost during the 
turn is a small number (~2x 10 -8 ) for oxygen burning, 
and is roughly the inverse of the Peclet number. Even 
within the narrow braking layer, there is little heat flow 
by radiative diffusion during oxygen burning. 

This discussion underestimates mixing because it ig¬ 
nores turbulent fluctuations (1 12.61) : larger fluctuations do 
more mixing than average, and mixing is irreversible. 
Turbulent kinetic energies fluctuate by factors ~ 2, 
so the mixing estimates should be increased accord¬ 
ingly. Flow velocities do not go to zero at the con¬ 
vective boundaries, but become small and oscillatory 
(iPressl IT oSlfc iPress fc Rvbickil 11981 1: iGoldrcich fc Kumarl 

ll99f)HGoldreich. Murray, fc Kumarll 1 0941) . As convective 

plumes hit the boundary, and rebound, the boundary 
moves in response; how elastic this is depends upon heat 
flow (the Peclet number). 

This “adiabatic” limit breaks down as the turnover 
time Ttum ~ b/u approaches the radiative diffusion time 
for the turn Tdiff ~ b 2 /Xc. For larger radiation mean- 
free-paths, the Peclet number decreases. No sharp tem¬ 
perature gradients can persist. This gives an “inelastic 
collision” of the flow with the boundary. This is the case 
for stars in photon-cooled stages of evolution; even with 
relatively large Peclet numbers for the whole convective 
region, the narrow boundary layers may still have sig¬ 
nificant energy flow by radiative diffusion. The previous 
discussion of the effect of excess pressure P' still holds, 
but because of thermal diffusion P' becomes increasingly 
dominated by density excess p' rather than the temper¬ 
ature excess T'. 

The red giant model of lViallet. et all (|2011l) provides 
an example of a boundary layer (the bott om) in_which 
there is significant radiative diffusion; iViallet, et all 
(120131) analyze this in detail (their § 4.6). As the bound¬ 
ary is approached from above, the down-flows are ac¬ 
celerated by pressure dilatation. These down-flows have 
an entropy deficit, so that they are heated by radiative 
diffusion from the surrounding material. In the braking 
region, compression causes a “hot spot” to develop. The 
flow is turned to a non-radial direction, and is now cooled 


by radiative diffusion (see Fig. 7 in lViallet. et al.lf2013f) . 

Such behavior differs from that obtained by present 
stellar evolution algorithms. The turning of the 
down-flow forces the mixed region to extend beyond 
that implied by the Schwarzschild criterion, and 
heating/cooling by radiative diffusion modifies the 
structure. While modest, such differences can be 
important for detailed models. In compensation for 
such changes, a standard solar model requires less 
opacity to have the same convection zone depth; this 
implies a lower metallicity. These changes in the solar 
model provide a means to reduce the disagreement 
with helioseismolo g v (IChrist ensen-Dalsgaar d. et al. I 
12011 iZhang et all 12012b iDcng fc Xional (l2008i) 


gave a justification for compositional s mooth¬ 
ing, as did sim ulations (IMeakin fc Arnettl I2007bt 
IViallet. et al.ll2013D. The thermal cha racteristics needed 
(IChristensen-Dalsgaard. et al. I 1201 If) follow from the 
analysis given above, which was not designed for the 
solar problem, and involved no solar or stellar calibra¬ 
tion. A more physically-correct convective boundary 
condition tends to improve agreement with abundances 
inferred from 3D stellar atmospheres (lAsplundl 120051) 
and the standard solar model. 

If heat flow processes are included, the “inelastic col¬ 
lision” with the boundary allows the loss of heat so that 
the entropy decreases for the downward flow, enhancing 
the downward acceleration. This effect tends to drive 
motion in convective envelopes. Heating at the bottom 
also tends to drive convective flow. However, c ooling at 
the bottom (as with URCA-shells. lArnettlll996h or heat¬ 
ing at the top (downwardly entrained, burning fuel) both 
tend to halt the flow. Such haltin g processes can cau se 
convective zones to split (jMocak. Siess fc Mullerl[20lTh . 

There may be observational evidence supporting this 
description of boundaries of convection which are deep 
in stellar interiors. Detection of g-mode pulsations in 
subdwarf B (sdB) stars allows an asteroseismic estima¬ 
tion of the size of the He-burning cores, which are sig¬ 
nificantly larger than predicted by the Schwarzschild 
criterion and standard stellar evolution theory (see 
iSchindlcr. Green fc Arnett] 120151 for discussion and ref¬ 
erences). Similar issues apparently are general for core 
heliu m burni ng stars observed by Kepler (iMosser et al.l 
l2014t iConstantino. Campbell fe Lattanzioll2014jh 

Finally, the origin (r = 0), in a ID stellar evolution¬ 
ary code using MLT, is a boundary as well. The use of 
Eq. [5] with adequate zoning implies that the convective 
velocity becomes very small due to symmetry (deriva¬ 
tives go to zero at the origin). This is a false braking 
layer caused by MLT being a local theory. Use of Eq. fal¬ 
lows flow through the origin provided that a counter flow 
gives conservation of linear momentum (e.g., a toroidal 
roll). At the origin in a turbulent convective core, this 
projects onto ID as a finite rms velocity, with a zero ra¬ 
dial gradient. MLT has problems with velocity at r = 0. 


4. SUMMARY 

We have brought more precision to the discussion of 
stellar convection by the use of 3D simulations of suf¬ 
ficient resolution to exhibit truly turbulent flow and 
boundary layers. The price paid is that we must re¬ 
place the unresolved turbulent cascade by Kolmogorov 
theory (ILES approximation), and the chaotic behavior 
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of an integral scale roll of Lorenz by a steady-state av¬ 
erage. We use RANS averaging to make 3D simulation 
data concise, and use 3D simulations to give RANS clo¬ 
sure. Solution of the RA NS equatio ns, using only the 
significant terms dMocak. et all 1201411 . is the full 321D 
procedure. 

This approach gives us a quantitative and precise foun¬ 
dation, based upon turbulent solutions of the equations 
of fluid dynamics. These numerical solutions have nu¬ 
merical limitations, which we have discussed. We find 
that the actual sub-grid dissipation in our simulations 
is automatically well approximated by the Kolmogorov 
four-fifths law. 

As a simpler first step, which addresses some of the 
worst errors of MLT, we focus on the acceleration equa¬ 
tion for the turbulent velocity. This makes the theory 
non-local, time dependent, and produces boundary lay¬ 
ers. It is almost identical to the equation developed 
from experimental study of the Rayleigh-Taylor insta¬ 
bility (RTI), indicating a close connection with plume 
models of convection; simulations also suggest this con¬ 
nection directly. Further development would entail use of 
RANS analysis to better deal with turbulent fluctuations 
f H2.6l and l2.7D . 

Even within the framework of the simple acceleration 
equation, there are several indications of how current 
practices in stellar evolution could be improved. The 
least drastic change involves diffusion: artificial diffu¬ 
sion 1< 13.2[) should be used with caution in situations in 
which real diffusion fi !3.3jl operates, because of distortion 
of the gradients which drive real diffusion (both artifi¬ 
cial and real diffusion have second-order spatial deriva¬ 
tives). The discussion in H3.8l gives a more realistic way to 
treat “overshooting”, and at the same time, removes the 
need for an imposed bou ndar y condition (Schwarzschild, 
Ledoux, or Richardson; H3.5D . The fluctuations in pres¬ 
sure discussed in < 13.81 will cause wave motion which will 
drive mixing in semi-convective regions on a dynamical 
timescale, far faste r th a n the thermal timescale con ven¬ 
tiona lly used (e.g., lLangcr. El Eid. fc Frickd (119851 ): see 

m- 

For use in stellar evolution this approach requires one 
more differential equation (for velocity, in addition to 
the traditional four, e.g., r, L , T, and p) and addi¬ 
tional coupling terms in the usual stellar evolution differ¬ 
ential equations (turbulent heating in the energy equa¬ 
tion, and ram pressure in the hydrostatic equation). The 
additional demand upon computational resources is not 
large. We use the convective flow velocity u and the 
super-adiabatic excess AV as separate variables, reflect- 
ing the fact that they ha ve different correlation lengths 
dMeakin fc Arnet^l2007blb We check that the simplified 
dynamic model does capture the numerical results of 3D 
as expressed in the RANS formulation. This approach 
is not calibrated to astronomical data, but predictive, 
being based on simulations and laboratory experiment. 
The simple 321D approach includes the Kolmogorov- 
Richardson turbulent cascade, and allows connections to 
past and future numerical simulations as a natural con¬ 
sequence. 

4.1. The future 

The enormous simplification, from 3D turbulent simu¬ 
lations requiring terabytes of storage down to a single 


additional ordinary differential equation (e.g., Eq. Ell- 
means that much is missing. For some applications 
the missing items may be important. One might use 
the RANS equations directly in a stellar evolutionary 
code , with 3D simulations to guide closure dMocak, et alj 
I2014D . We have presented a step toward that goal. Alter¬ 
natively, one might add to the simple 321D as needed, 
using new models guided by RANS results. Probably 
both paths should be followed, given the complexity of 
the problem. 


4.1.1. 321D algorithms 

We have refrained from offering detailed algorithms be¬ 
cause we believe that there may be a variety of useful 
ones, tailored for existing stellar evolution codes, and to 
be modified by developing insight. This is not a finished 
subject. A skeleton algorithm should include: 


1. velocity from an acceleration equation (Eq.[6l H2.3D . 

2 . boundary physics: turning, damping, mixing and 
shear ( ^3.81) . 

3. fluxes of enthalpy and composition f l|2.4l and 1 13.81) . 


4. non-locality in velocity: turbulent kinetic energy 
flux and ram pressure 1 12.31) . and 

5. turbulent heating of background by Kolmogorov 
cascade (Eq. Q]). 


Our first priority is to implement these ideas in 


in TYCHO (Liebert. et al. 120' 

11 ), and plan to mi- 

grate to MESA (|Paxton. et al. 

:20nl (20131). MON- 

STAR flCampbell & Lattanzio 
20101). CENEO Hones, et al. 

20081: Dohertv. et al. 
20151). and FR.ANEC 


plementations in other codes. 


4.1.2. Further simulations 

New simulations to better quantify the bound¬ 
ary physics are in progress (|Campbell. et al. I 120151 
iCristini. et al. I l~2015l) . This approach, unlike MLT, is 
generalizable i n principle to include ro tation and MHD 
(|Maederlll999t iMaeder fc Mevnetll2000D because it starts 
with full 3D equations. For example, rotationa l terms 
are im plicit in t he vector form of Eq . [6] see also iBalbusI 
(|2009l ) : iFeathcrston fc Mieschl (120151 ). 

4.2. Implications 

Because of the fundamental importance of convection 
in stellar evolution theory, a replacement for MLT will 
have implications for many areas throughout astronomy 
and astrophysics. A few of the most striking are: 

4.2.1. Helioseismology 

Convective boundaries with low Peclet number will 
be smoother, which reduces the disagreement be¬ 
tween helioseismology and solar model predictions; see 
iChristensen-Dalsgaard. et al. 1 120111 : iZhang et al. 1 120121 
and 

The corrected boundary conditions for convection 
will place the composition gradient further beyond the 




















































Beyond MLT 


17 


Schwarzschild zero condition PI3.8D . requiring a lower 
opacity below the mixing boundary to get an acceptable 
solar model. This may be attained by a lower metallicity, 
which will reduce the disagreement between solar models, 
and solar abun dances determined from 3D atmospheres 
( lAsplundl 1200511 . The combination of these two correc¬ 
tions will shift the standard solar model problem toward 
the Asplund abundances. 

4.2.2. Asteroseismology 

These modifications beyond MLT bear on many dis¬ 
crepancies between asteroseismology and stellar evo¬ 
lution theory. Some examples: application of bet¬ 
ter convective boundary physics will produce larger 
He burning cores in sdB stars, and reduce the large 
discrepancy between the asteroseismolo gy determ ina¬ 
tion o f core sizes and ste l lar m o dels (ICharpinet et all 
19971: IVan Grootel. et al. I 201C iBloemen et al.l 12014 
Schindler. Green fe Arnettl 12015 1. Similar issues ap¬ 
parently are general for core he lium burning stars ob¬ 
served by Kepler dMosser et al.1 12014H . The discrep¬ 
ancy i n mixed modes in normal CHeB (“red clump”! 
stars (iBildstens. et al. I 120121: iMontalban. et al. I 12013c 

Stello. et al. 1 120131 iConstantino. Campbell fc Lattanziol 
2014fl will be affected. 

4.2.3. Convective boundaries, nucleosynthesis yields and 
pre-supernovae 

The nature of convective boundaries is affected by ra¬ 
diative diffusion, so that they differ for neutrino-cooled 
stages of nuclear burning. Calibration of convection for 
late stages, from stages dominated by photon-cooling, 
requires re-evaluation. Detailed estimates of stellar nu¬ 
cleosynthesis and stellar struc ture based u pon an algo¬ 
rith mic diffusi on scenario (e.g.. IWooslev fc WeavefllT995l : 
iWooslev. Heger. fc \Veaveill2002l) are not confirmed, and 

require re-examination. 

Whil e the general features of nucleosynthesis yields are 
robust (lArnettlll996ll . detailed abundances depend upon 
details of mixing and convection. Nucleosynthesis from 
lower mass stars is also affected: asymptotic giant branch 
(AGB) stars do not have a third dredge up without “over¬ 
shoot”, which is a convective boundary problem. This 
dredge up is crucial for s-process nucleosynthe sis (it pro¬ 
vides a neutron source, lLattanzio. et al.ll2014T) . 

Driven by neutrino cooling, nuclear burning in stars 
prior to core collapse is vigorous, and in turn drives vig¬ 
orous convection. Convective velocities increase as evo¬ 
lution proceeds. The nuclear energy generation is, on 
average, in balance with the turbulent dissipation at the 
Kolmogorov scale, so e nuc ~ u 3 /£, which relates the nu¬ 
clear energy generation rate, the average convective ve¬ 
locity, and the depth of the convective zone. Velocity 
fluctuations are large dMeakin fc Arnettll2007bfl . Super¬ 
nova progenitor models which are ID can represent av¬ 
erage properties, such as convective speed, but not the 


amplitude and phase of the (large) fluctuations of those 
properties. Re alistic progenitor model s should be dy¬ 
namic and 3D (|Arnett fc Meakinll2011allbl ') if they are to 
be used for accurate core collapse simulations. 

4.2.4. Core collapse 

The size and structure of progenitor cores affects the 
possibil ity of producing explosions in cor e collapse simu¬ 
lations dCouch fc Qttll2013tlArnettll2014[) . The predicted 
size and structure of such cores depends upon the physics 
of convection used in the stellar evolution codes. Detailed 
scenarios for pre-super nova structure, c oll apse and explo¬ 
sion, such as found in IWooslev. Heger. fc WeaveJ (120021 ) 
for example, are not robust, and may require revision 
when better treatments of mixing are applied. The va¬ 
lidity of calibrating neutrino cooled convection on pho¬ 
ton cooled stages of evolution is questionable due to the 
large difference in Peclet number. Even the size of the He 
core is uncertain with present algorithms (lLangcrlll991l 
i2oHh and will be affected by better treatment of con¬ 
vection and convective boundaries. The theoretical ap¬ 
proach to turbulence used ab ove can also be applied t o 
the core collapse process itself dMurnhv fc Meakinll2011[l , 
giving insight even for 3D simulations which are presently 
under-resolved due to computational limitations. 
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APPENDIX 

THE CONVECTION EQUATIONS 

We develop the fluid equations in an inertial frame (lLandau fc Lifshi~t3ll959f) . We begin with a general formulation, 
and transition to a specifically spherical (r, 9, <f) choice of coordinates for application to stars. We will decompose 
variables into a background part and a fluctuating part, e.g., for pressure P = Pq + P'. Our procedure is chosen for 
stars in which the background is hydrostatic and spherically symmetric, so that VPq = — pog = —g/Vo. 
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Baryon Conservation 

The vector form of the continuity equation (lLandau fc Lifshit3ll959f l is 

dp/dt + V-pu = 0, (Al) 

where p is the mass density and u is the fluid velocity. In the incompressible limit, for a steady flow, the net flux 
of mass into a region equals the mass flux out. In thin boundary layer, perpendicular to the radial direction r, the 
average velocities must satisfy 

du r /dr = —2 duh/dh, (A2) 

where h is either of the symmetric transverse coordinates (i.e., locally cartesian), to avoid changing the density (as 
s een i n the E ule rian fr ame). 

IViallet. et all (2013) show (their Eq. 28), that for fluctuations against a steady background, 


Vu'= 


H n 


(A3) 


where H p is the density scale height, and u' r is the radial component of the velocity fluctuation. This approaches 
zero (the incompressible limit) for shallow, subsonic convection (large density scale height and small radial velocity 
mach number, u' r -C s, where s is the sound speed). This velocity “dilatation” is due to the vertical motion in 
the background stra tification and becomes an important component in convective driving i n deep convection zones 
(IViallet. et al.ll2013f). N otice that rising plumes ( u r > 0) expand and falling plumes contract dStein fc Nordlundlll989l : 
iMeakin fc Arnettll2010D . 


Momentum Conservation 

The vector acceleration equation (Eq. [6]) is 


du/dt + (u ■ V)u = B — u/r 


(A4) 


where u is the velocity, r 

If 


\u\/£d with £d is the Kolmogorov damping length, and the variable B is defined as in 112.31 


B = --S7P- g, 
P 


(A5) 


where P is pressure and g is gravitational acceleration, then Eq. IA4I is a Navier-Stoke s desc ription of the largest 
scales of turbulence, with a simplified damping term which is consistent with iKolmogorovI (Il962f) . Note that the usual 
formulation of hydrostatic equilibrium in stellar evolution theory is some variant of the condition B = 0. Projecting 
Eq.ESI onto the radial coordinate, we have 


du r /dt + u r du r /dr = ——dP/dr — g — u r /r. (A6) 

P 

Th e full equations in spherical coordinates are shown in S15. lLandau fc Lifshitd (|1959f) (see also lMihalas fc Mihalasl 
119841 for a detailed discussion), with the bare viscosity terms rather than Komogorov’s expression for integration of 
the turbulent cascade. In tensor form the momentum equation is 


dui/dt + Ukdui/d Xk = — dP/dxi — gi 

P 

1 d r /dui_ duk_ _ 2^ dui \i d / ,du L \ 
p dxk V dxk dxi 3 lk dxi ). dxi V dxi) 


(A7) 


Kolmogorov’s four-fifths law ( Frischl[l995l) stat es an amazing simplification, that integration over the turbulent cascade 
reduces the last term in Eq. [AT] to —u/r ('Eq. IA4I) on average, ignoring boundary effects (see ©. 

To illustrate how turning happens at boundaries, it is sufficient to consider the simpler case of flows with 6 and 
4> length scales small compared to r, so the transverse dimensions are quasi-cartesian (the inertial terms in 1 jr are 
neglected; for convective cores, the more cumbersome full equations are needed because r cannot be large near the 
origin). Then the two transverse components are symmetric in this approximation and satisfy 


duh/dt + Uhduh/dh = —~dP/dh — Uh/r, 

P 


(AS) 


where dh is rdd or r sin 6d(f>. We consider finite fluctuations about a static backgroun d, so that we substitute p = po + p 1 
and P = Pq+P'- We ignore variations in g (the Cowling approximation. ICoxHl980H . Using —dPo/dr = p^g, the radial 
equation becomes 

du r /dt + u r du r /dr = — f———- 

\p 0 +p'J Po + P 


-dP'/dr — Ut/t. 


(A9) 
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Convection is often described using only the buoyancy term; the pressure fluctuations are taken to be small, of order the 
mach number squared. However, near boundaries the pressure fluctuations pr ovide the tangential acceleration which 
is necessary to turn the flow, and should not be neglected (see lNordlundlll985l) . The buoyancy term acts through the 
density fluctuation p', and only in the direction parallel to the gravity vector. The transverse equation is 

duh/dt + uhduh/dh = --— -dP'/dh — Uh/r. (A10) 

Po + P 

Note that the radial and transverse equations are coupled primarily by the pressure fluctuation term P\ but also by 
u/t, because r = i d / |w| where \u \ 2 = u 2 = u 2 + 2 u 2 (turbulence damps regardless of orientation of the large scale 
flow). The fluctuating pressure near convective boundaries insures the generation of waves. 


Energy Conservation 

Following lLandau fe Lifshitzl (11959H . §6, the equation of energy conservation is 


d_ 

dt 


1 


- pu 2 + pE + 


-V- 


pu{-u 2 


■W + <j>) 


+ T 


dpS 

~dT' 


(All) 


where <f> is the gravitational pot ential and g = — V d>. If taken to both the steady state and adiabatic limits, this 
becomes the Bernoulli equation (lLandau fe Lifshitzlll959h . The entropy change equation may be written as 


^PS 

T nr = '*■ 


F r 


(A12) 


where e nuc is the net heating from nuclear and neutrino reactions, pt V i SC is the Navier-Stokes viscous heating term as 
modified by Kolmogorov’s four-fifth’s law (see Ea.[Tl lA4l and lA7l) . and F rad is the energy flux due to radiative diffusion. 
The viscous term is missing from MLT and the Euler equation. Most of the turbulent kinetic energy resides in the 
largest (integral) scale, while turbulent heating occurs at the small (Kolmogorov) scale. Then tturb = u • u|it| jld is 
the Kolm ogorov heating from the turbulent cascade, and TdpS/dt , pe nuc and F ra d are now the appr opriate RANS 
averag es (IViallet. et al.N20l3lh One req uirem ent for Berno ulli’s equation to be vali d, as assumed in iPasetto. et al.l 
(12014!) (see H2.8|) . is that the RHS of Eq. IA12lmust be zero (lLanda u fc Lifshitz 119591 Ch. I). This is found not to be 


gener ally true, either in the 3D s imulations dViallet, et alJl2013t iMocak, et al.l 2014*1) . or experimentally in turbulent 
flows dTennekes fc Lumlevl Il972t iDavidsonl 120041) . Heating is an essential feature of 3D turbulence, which converts 
large scale, ordered velocities to disordered ones. 
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